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SYNOPSIS 

A. KRISHNAN 
Ph.D. 

Department of Electrical Engineering 
Indian Institute of Technology , Kanpur. 

November, 1978 

STUDIES IN NONLINEAR DISCRETE TIME SYSTEMS 

For nonlinear systems, described by either nonlinear 
differential or nonlinear difference equations there is unfor- 
tunately no general procedure for obtaining solutions m closed 
form. In fact, when a physical problem leads to nonlinear 
differential (or difference) equations, one is most often 
content with a qualitative description of the solution together 
with a numerical approximation. Most engineering situations 
are inherently nonlinear, and can be successfully linearised 
only for particular operating condition in order to apply 
known linear techniques. However, a considerable number of 
useful engineering problems can not be adequately linearised. 

Both weakly and strongly nonlinear continuous time systems 
have been investigated by various persons. However, little 
work has been reported on the analysis of discrete time systems, 
in particular nonlinear discrete systems. Since I960 discrete 

i 

systems have gone through a rapid change m the control area 
and the principle of optimal control, the concept of state 
variables, stability studies by Laapunov's method etc. have 
been used for discrete systems study. The well known technique 



of Krylov - Bogoliubov averaging method has also been applied tc 
a restricted class of Nonlinear Discrete Time Systems to obtain 
approximate solution. 

i 

In this thesis, nonlinear difference equations are consi- 
dered and a detailed analysis is carried out for both autono- 
mous and nonautonomous situations. Nonlinear difference 
equations arise commonly in the analysis of various discrete 
time problems such as sampled data control systems, digital 
filters, biological oscillations etc. Analysis of nonlinear 
difference equations m the time domain have however been 
studied to a much lesser degree than their counterpart - non- 
linear differential equations. The major contributions m this 
thesis are listed below : 

1. A new method, known as discrete multiple scale pertur- 
bational technique, for analysis of nonlinear difference equa- 
tions is introduced, which is sufficiently general to consider 
both autonomous and nonautonomous descriptions. This method 
is based on the introduction of two separate independent time 
scales - a fast time scale and a slow time scale. This 
technique has been used extensively to study a class of 
differential equations (both ordinary and partial), but there 
has been no instance of its application to discrete time system 
oxcept the paper by Hoppensteadt et al. dealing with a system 
of difference equations in matrix form, m which the discrete 


time model is obtained through, the well known Taylor series 
expansion. The proposed multiple scale perturbational tech- 
nique, on the other hand is developed using the known proper- 
ties of finite difference operators. 

2. The discrete multiple scale perturbational technique has 
been applied to a class of nonlinear difference equations to 
obtain both qualitative and quantitative information, linear 
and nonlinear difference equations under free and forced 
situations are analysed using the proposed technique. Super/sul 
harmonic oscillations under strongly forced condition are also 
investigated through the proposed technique. 

3. iJ -'he multiple time method is applicable only for a polynomia. 
kind of nonlinearity . The remaining part of the thesis is 
devoted to the analysis of nonlinear dil Lerence equations 
with saturation type of nonlmearities which are common in 
digital filters. These nonlmearities are due to adder over- 
flow and quantization during addition or multiplication of two 
or more signals. The present study is limited to second order 
digital filters (which are in fact basic blocks for higher 
order systems) with overflow saturation nonlinearity . The 
present study has focussed on the following new investigations 

(i) the existence of limit cycles with dilforent periods for 
parameter values outsid e the stability triangle. 

(li) the location of a region outside the stability triangle in 
the parameter plane in which the output sequence dies down 
to zero monotomcally . 
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(111) under forced situation, conditions arc derived for the 
existence of jump phenomenon as well as nonlinear sus- 
tained oscillations when the paramo ter values, the magnituc 
and the frequency of the input function are known before- 
hand. 

a chapterwisc summary of the work reported m this thesis 
is given now ; 

The first chapter provides an mt reduction to the non- 
linear systems and nonlinear discrete lime systems m particular. 
A survey of some of the earlier work timing to this thesis 
10 presented. A brief description of s ecu, id order digital 
filters, the nonlinear phenomenon mtrod icoO auo to finite 
wordlength registers and their effect on normal operation 
of the filter are discussed. The scope .010 objectives of the 
thesis arc also outlined. 

The second chapter is concerned wiou nonlinear difference 
equations with weak polynomial typo of nonlinearity under free 
and weakly forced situations. The proposed discrete multiple 
scale perturbational scheme is applied to linear as well as a 
class of nonlinear difference equations and the deduced results 
are compared with exact solution obtained by computer simulation, 
for weakly nonlinear systems, a Vander pol typo model is con- 
sidered and its limit cycle oscillations arc investigated. The 
stability property is studied by a variational approach. The 
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discrete version of the Duffing equation under free and 
weakly forced situations are also taken as illustrative 
examples and the response characteristics such as variation 
of amplitude with input frequency are obtained. For weakly 
forced and weakly damped case, the concepts like vertical and 
horizontal tangents m the steady state response characteri- 
stics are also studied. 

Chapter three deals with the important concept of super/ 
sub harmonic oscillations m nonlinear discrete time systems 
under strongly forced situation and the analysis is carried 
out through the multiple scale per turbat Lonal technique as 
well as by harmonic balancing method. A Duffing type 
equation is considered and the possible cup or /sub harmonic 
oscillations are obtained. 

In the fourth chapter a saturation type of nonlinearity 
is considered which is common m digital fillers. For a second 
order digital filter with overflow saturation nonlinearity 
under force free condition analytical c.ipi’ocsions are derived 
to obtain different regions in the a~b pirunctcr plane for 
sustenance of different periods of limit cycle oscillations. 
Besides the stability investigation anot i „r interesting 
contribution is to locate the region outside the stability 
triangle m the a-b parameter plane m w-'ich the output 
sequence decays to the zero solution mono conically for a 



particular set of initial values of the filter variable. 

The region in the parameter plane and tne initial condition 
loci are obtained for such a monotonic re op one c . 

The fifth chapter focusses on the investigation m digital 
filter under forced situation, known as ' jump phenomenon' . 

Jumps exist m second order digital filters, for a specific 
filter parameter values and the magnitude of periodic input 
function, whenever there is a disturbance in the filter variable 
For a given input frequency , conditions are derived in tens 
of tne filter coefficients and the Magnitude of the input fun- 
ction for the existence of a jump or nonlinear sustained osci- 
llation. Location of regions in the a-b parameter plane for 
existence of such nonlinear phenomena are also shown by graphi- 
cal construction of the derived conditions. In addition to 
the above nonlinear phenomena m forced digital filter, the 
possibility of the existence of subharmonic oscillations are 
also investigated. 

The work reported and the contributions made m this 
thesis arc reviewed m the sixth chapter. This chapter 
concludes with an assessment of the scope for further research 


work in this area 



CHAPT3B. 1 
INTRODUCTION 


1*1, General Discussion : 

The dynamic analysis of a physical system is usually a 
description, in mathematical form, of the physical properties 
of the system. This form is often called tho mathematical 
model of the system. Most mathematical models of physical 
systems are integro-diff ercntial equations or difference 
equations, or a combination of those equations. When these 
equations are linear the system is said to be linear, which 
are in turn classified into those with constant and variable 
parameters. In linear systems with variable parameters, the 
parameters are functions of the independent variable, genera- 
lly the time. These systems are also known as parametric 
systems, a system is nonlinear if its behaviour can not be 
described by linear equations. Physically this means that 
the parameters are functions of dependent variable and its 
derivatives. The most significant difference between linear 
and nonlinear equations, from an engineering point of view, 
is that linear equations obey the laws of superposition, -a* ^ h ,k> 
while nonlinear equations do not. 

Nonlinear equations, which describe the behaviour of 


nonlinear systems, cannot as a rule be solved explicitly, 
unlike linear systems with constant coefficients. Among 
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nonlinear systems the separate class of self-excited systems 
have been the subject matter of extensive research. This 
class includes systems that can produce oscillations without 
an external force. In general nonlinear systems are an 
indispensable part of most engineering applications, and their 
study is of much importance. Various techniques for inves- 
tigating the solutions are used to aid in their study. These 
techniques are numerous and a brief survey is given in the 
sections 1.2 and 1.4. The next two sections describe some 
salient features of nonlinear continuous time and nonlinear 
discrete time systems. 

1.2 Honlinear Continuous Time Systems ; 

Nonlineanties exist in almost all physical problems . 

For instance many problems in aerodynamics and hydrodynamics 
involve nonlinear equations |_1] . In electronics a common 
nonlinear device is the triode oscillator. Appel ton and 
Van der Pol |_2,3J first derived a mathematical model for 
this circuit, des cribed by nonlinear differential equation. la 
Van dcr P^l obtained solutions for this equation by 
the isocline method. In magnetic circuits, saturation of 
the iron core is a nonlinear phenomenon. In mechanical 
systems, the mass on a nonlinear spring is a nonlinear 
system which was studied by Duffing and others in early 1920. 

In general a nonlinear system can be analysed by 
(a) qualitative methods and (b) quantitative methods. In 
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qualitative analysis, the complete information can be obtained 
without having to solve the given equations. The phase piano 
analysis and describing function techniques are some of the 
techniques employed. The quantitative methods are analyti- 
cal in nature. The study of nonlinear systems through these 
methods has been a subject of much study since the time of 
Poincare. The analytical techniques used m nonlinear systems 
analysis have been developed as a result of work done m 
nonlinear mechanics and nonlinear circuit analysis. However, 
at present it appears that no generalised theory of solution 
can be formulated for all nonlinear systems. 

The behaviour of a nonlinear system is quite different 
from that of a linear system. Unlike linear systems, a 
nonlinear system may be stable for small inputs and may be- 
come unstable if the input exceeds a certain level. The 
nonlinear system may respond in certain other ways which 
are peculiar only to this type of system. Limit cycle 
oscillations, super /subharmonic oscillations, froquency 
entrainments and jump phenomena are some of the examples 
for the peculiar behaviour of nonlinear systems. These are 
briefly explained as follows : 

1.2.1 Limit cycle oscillation : 

This form of oscillation in a nonlinear system is a 
steady state oscillation, usually periodic but not necessaril; 
sinusoidal. In the phase plane the limit cycle oscillations 
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have closed curves known as periodic mjccancs. Stability 
of such periodic trajectories arc studied by considering a 
small variation on either sides of the closed trajectory. 
Periodic trajectories which arc asymptotically stable from 
both sides are called stable limt cycles. Those which are 
unstable from both sidos arc called unstable limit cycles. 

The periodic trajectories which arc stable from one side and 
unstable from the other side arc known as corn is table limit 
cycles. Some nonlinear systems may have multiple limit cycles 
whose trajectories arc closed isolited curves. 

1.2. 2 Super and subharmonic oscillations : 

In a linear system where input is a sinusoidal forcing 
function with frequency w , the steady state output will also 
oscillate at frequency w , except in the special case whore 
the system has two or more repeated poles on the imaginary 
axis. In some nonlinear systems , on the other hand* for the 
same forcing function, the output will have a steady state 
oscillation that is comprised of the fundamental component 
at frequency u> along with components at frequency pu /q 
whore p and q are positive integers . When the frequency 
component p w /q is predominant then tlio steady state oscilla- 
tions are classified as follows : 

(l) Super or ultraharmonic oscillations,* if q = 1 and p- 2, 
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(li) Subharmonic oscillations,' if p = 1 a^d q 2, and 
(ni) Ultrasubharmonic oscillations,* if p and q are relative 
prime integers f (p and q are said to be relative prime 
integers if their greatest common divisor is unity.") 

1.2.3 Frequency entrainment ; 

If a nonlinear system has a natural limit cycle osci- 
llation of frequency and if it is forced by an external 
input at frequency U2 » caay be expected that components 
at both frequencies will appear at the output. In some 
cases , however if w 2 ^- s c l° se to w ,, the output component 
at will vanish leaving an oscillation at &>2 only* This 
phenomenon is known as frequency entrainment. The range of 
values for corresponding to to ^ for which this phenomenon 
occurs is called the zone of synchronization. When w 2 is not 
close to w -j.* both components will generally be present. A 
considerable amount of literature concerning the frequency 
entrainment £3 available which includes tlio text book 
by Hayashi [4]. In this book the possible occurrence of 
highcr/subharmonic entrainments are also discussed. 

1.2.4 Jump phenomenon ; 

Yet another phenomenon oceunng m some nonlinear 
systems is discontinuous jumps m output magnitude. This 
happens for smooth and continuous variation of either fre- 
quency or amplitude of the forcing function. For either 
*wtiere their greatest common divisor is unity. 
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variation, the output jumps from one value to another with 
widely differing magnitude. This has been analysed m great 
detail m many text books including Stoker (_5], Gunningham [6] 
and Ku [7]. 

The existence of these and other phenomena that are not 
found m linear systems makes the field of nonlinear systems 
quite complicated for study. The theory of nonlinear systems 
has not been developed to the point where one can say very 
much about the solution the system from the equation it- 
self, as is possible with linear constant coefficient equa- 
tions. Therefore, each nonlinear equation pres exits a new 
problem. Now a brief review of the met nods available to 
analyse the nonlinear systems is presented below. 

The methods for finding solutions l'or nonlinear differen- 
tial equations or differential equations with varying coeffi- 
cients t generally are more diLXicult and the solutions 
less satisfying on account of the degree of approximation 
involved, than those for linear coxistant coefficient differen- 
tial equations. Only m few cases can exact solution in 
terms of known functions be found. Usually, an approximate 
solution is possible, and this solution nay be valid with 
reasonable accuracy only for certain range of values of the 
independent variable. When the equations are nonlinear, 
solutions by analytical methods are possible only if the 
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nonlinearity is small and even then available solutions are 
generally approximate, The system described by a small 
nonlinearity is generally known as weakly nonlinear systems, 
Equations with a considerable degree of nonlinearity are 
known as strongly nonlinear systems, These systems can be 
solved generally by numerical or graphical methods. 

Several procedures have been discussed by Ounnigham [6], 

Mela chlan [8], Andronov et al [9], Gibson [10] and many 
others for obtaining the approximate solution for a nonlinear 
differential equation by analytical approach. The method of 
slowly varying amplitudes* power series methods and pertur- 
bational methods are commonly used to analyse suoh systems. 

The method of slowly varying amplitude© is widely used in 
the theory of oscillations. This method was first introduced by 
Van dor Pol [11], who had considered a number of problems 
on transient phenomena in valve oscillators and other systems. 
This method has been further extend od by iDcylov and Bogoliubov 
[12] among others to study nonlinear systems* The technique 
developed by Khylov-Bogoliubov has been amplified and justi- 
fied by Bogoli^ov and fiLtropolsfci [13] who have developed 
the powerful technique known as itrylov-Bogoliubov-Mitropolski 
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is sufficiently snail, a first corrected solution by 
applying either power series or perturbation methods may 
give sufficient accuracy. If the degree of nonlinearity is 
large, It is sometimes possible to obtain better accuracy 
by applying those methods repeatedly* Further repeated 
application of these methods are possible in theory, but in 
practice the mathematical computations become comp li cat ad. 

The perturbation technique is applicable to equations with 
a small parameter associated with nonlinear terms, kn appro- 
ximate solution is found as a power series, with terms invol- 
ving the small parameter raised to successively higher powers. 
If the magnitude of the parameter is small enough, the first 
few terms of the series are sufficient to give a fairly 
accurate solution. 

Most of the available analytical methods are limited 
in their application to the weakly nonlinear systems. The 
perturbation technique is most widely used in the field of 
nonlinear systems, fhis method accounts for the effect of 
nonlinearity by additive correction terms and is also rela- 
tively straiglitforward at least where the linear terms as 
well ae nonlinear terms appear in the differential equations , 
Where the solution is oscillatory, however, secular terms 
with indefinitely increasing magnitude may arise, and the 
procedure must bo modified to eliminate such term© [6]. 
fhis method my be applied successively in order to obtain 
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additional terms in the series solution and thereby to 
achieve better accuracy. Each successive application, how- 
ever, generally becomes tedious. 

When the solution of nonlinear equation is an oscilla- 
tion having its amplitude or phase changing with time, but 
doing so relatively slowly, the multiple tine perturbation 
technique is useful. 

The two time scaling which is a special case of multiple 
time perturbation technique has been used successfully for 
analysis of a class of second order differential and partial 
differential equations. The formulation of two time expan- 
sion scheme was initially given by Cole and Kevorkian [14] 
and they obtained uniformly valid asymptotic approximations 
for certain nonlinear differential equations. Kevorkian [15] , 
Cole [16], Carrier and Pearson [17] and Nayfeh [18] have 
discussed in great detail, the theory and the application of 
two variable expansion in analysing a class of nonlinear 
second order differential equations as well as partial diffe- 
rential equations. Rasmusson [19 J obtained uniformly valid 
zeroth order approximation for fchc general Lienard type 
equation with small damping by means of multiple time scaling. 
The two time scaling method has been extended to analyse a 
set of coupled as well as higher order systems. Yen and 
Kronauer [20] studied the response of three oscillators 
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coupled by a weeds nonlinearity through, this approach. Re- 
cently Tiwari and Subramanian [ 21] have successfully applied 
the two time scaling technique for third order nonlinear 
differential equations of various classes with polynomial 
kind of nonlinearities. Tiwari et al [22] and Krona uen- et 
al [25] have considered the application of two variable 
expansion scheme for determining super, sub and ultrasub- 
harmonic resonances in a class of nonlinear second order 
differential equations with multiple and single input condi- 
tions respectively, lick [24] and Jordan [25] have also 
considered the two time scaling approach to analyse the 
singular perturbation problems and Mathiou equations respe- 
ctively. In short, the nonlinear oscillatory systems have 
been analysed deeply by various persons using this techniques. 

During the past few years transformation techniques 
have also found considerable application in tho study of non- 
linear differential equations. Dasarathy and Srinivasan [26] 
used transformation methods to reduce a restricted class of 
nonlinear systems into equivalent linear systems. Jones 
and Ames [27] applied similarity methods to reduce the order 
of nonlinear differential equations. Musa and Srivastava 
[28] employed Ido's method of infinitesimal bra ns formations 
to generate some general classes of second order nonlinear 
nonautonomous differential equations which can be reduced 
to first order systems. This method is most general and 
applicable to higher order systems too. 
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In addition to the above approximate methods, a method 
of analysing nonlinear systems known as functional analysis 
approach has also been established in the literature. The 
functional analysis viewpoint provides a broad and powerful 
means for treating many problems m nonlinear systems as 
well as in other fields like filtering and optimal control. 
The book by Porter [29] gives a detailed introduction to 
linear functional analysis with engineering motivation. A 
brief introduction to functional analysis approach and its 
application to nonlinear systems has been given by Holtzman 
[30]. In [30] the functional analysis approach has been 
successfully applied to study nonlinear systems especially 
periodic solution of weakly nonlinear systems and stability 
analysis of nonlinear systems. 

1.3 Nonlinear Discrete Time Systems : 

Unlike the continuous time systems, in discrete time 
systems the dependent variable assumes values only at dis- 
crete time instants. Difference equations bear the same 
relation to the discrete time systems as differential equa- 
tions do to continuous time systems. As m the case of 
continuous time systems, the discrete time systems can also 
be classified as linear and nonlinear systems. The system 
is said to be linear if the mathematical model namely the 
difference equation is linear, that is the dependent variable 
and its differences appear only to the first power. If this 
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power is other than the first, or if the vanableepppear as 
product with one another, the equations are nonlinear. In 
general a difference equation is defined as follows. Let 
b denote a set of real numbers x. The set b may include all 
real numbers, or it may be restricted to integers, or non- 
negative integers, or only positive integers. Then if a 

function y(x) is defined over the set S, an equation relating 

2 

the values of y and one or more of its differences A y, A y, 

A u y is called a ’Difference Equation’ . The symbol & 

denotes the forward difference operator. The difference 
equation can also be represent ed in terms of central and 
backward difference operators [31,32,33 J. A detailed des- 
cription of these operators is given in Chapter two. 

Difference equations also occur in combined form with 
differential equations, particularly in the field of stochas- 
tic processes. This form is known as difference-differential 
equations. In difference equations the independent variable 
may be discrete (usually a nonnegative interger valued variabli 
or continuous. In the latter case the equation is sometimes 
called a functional difference equation. 

Difference equations not only play a role m their own 
right as direct mathematical models of physical phenomena but 
also provide the field of numerical analysis with powerful 
tools. By approximating differential equations with difforenc 
equations and thereby obtaining numerical solutions, one can 
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study approximate numerical solutions of nonlinear differen- 
tial and integral equations for which there are no known 
closed forms. 

Difference equations arise commonly 11 . analysis of 
various discrete time problems such as digital filters, 
sampled data control systems, biological oscillations etc. 

A survey paper by Jury and Tsypkin [34] indicates the spect- 
rum of research m discrete time systems over nearly two 
decades since 1950. Various applications are also cited m 
the above reference. 

In many fields of Engineering and Science the behaviour 
of certain systems may be modelled in terms of difference 
equations instead of differential equations. For instance 
Hsu and Yee [35 J have studied the behaviour of dynamical 
systems governed by a simple nonlinear difference equations. 
Lorenz [36] in connection with meteorology and fluid flow 
stability problems, May [37] for population dynamics problems 
and Hsu and Cheng [58] for dynamical systems governed by 
periodic differential equations, also obtained difference 
equations for their models. Thus the development of theory 
of dynamical systems described by difference equations 
seems to hold promise of having applications m many areas 
of research. 

As mentioned earlier, difference equations also exist 
in digital filters. Due to limitation in storing the signals 
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in registers, under certain circumstances a digital filter 
exhibits nonlinear operations. The nonlinearity introduced 
due to nonlinear operations is of discontinuous type and they 
can not be adequately represented by polynomials. Hence 
some special methods are to be adopted to analyse such systems. 
The following section is concerned with digital filters and 
the nonlinear phenomena associated with them. 

1.3.1 Digital filters : 

The nonlinearity associated with the discrete time 
system may be of a polynomial typo or piecewise linear e.g. 
saturation type. The latter type of nonlinearity is common 
in physical systems such as relay servomechanisms and digital 
filters. Digital filtering is one of the important tools 
for modern digital signal processing and it has found many 
applications in an increasing number of fields m science 
and engineering. A digital filter may be defined [ 39 ] as a 
computational process or algorithm which converts one sequence 
of numbers representing an input signal into another sequence 
of numbers representing an output signal, and m which the 
conversion changes the character of the signal m some pres- 
cribed fashion. 

Digital filters are discrete systems, with the addi- 
tional property that the signals at any instant can only 
have a finite number of values, usually represented in binary 
form. This makes it possible to store the signals in 
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registers, necessarily having a finite length. In spite of 
many advantages offered by digital filters there is an 
inherent limitation on the accuracy of these filters due to 
this finite wordlength. A survey and discussion of the effect 
of finite wordlength on the accuracy of digital filters is 
given by Liu [40] and Qppenhem et al [41] . In such systems 
the finite wordlength available for the signals affects the 
perfect operation and in most cases introduces nonlinear 
phenomena. 

The nonlinear phenomena introduced m digital filters 
are due to arithmetic operations, such as multiplications 
and additions performed on the signals. After these opera- 
tions, the wordlength will, in general, have increased. Assu- 
ming a fixed point representation [42, 43] a multiplication of 
an N bit word with another £ bit word results in a word having 
(N+IC) bits. Similarly addition of two words of N bits will 
result m a word requiring (N+l) bits for its representation. 
This increasing word length could be handled by providing lar- 
ger registers, which in most systems, would lead to very im- 
practical size. To avoid this kind of practical difficulty, 
often a wordlength reduction is applied. Generally two types 
of wordlength reduction namely quantization and overflow are 
applied. 

(l) Quantization ; 

This is obtained by either round ing a digital word 
if say N+K bits to N bits or supressmg several of the least 
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significant bits (magnitude truncation). This effectively 
results in a nonlinear digital filter. 

(ii) Overflow ; 

Overflow can occur after a multiplication or an addition 
and this affects the most significant bits. Consequently 
the overflow causes more severe signal distortion than the 
quantization. 

Both of these operations result in a deviation of actual 
output from the desired output. Because of these nonlinear 
operations limit cycle oscillations can occur under zero 
input situation and nonlinear sustained oscillations, jump 
phenomena, super /subharmonic oscillations are also possible 
under forced input situation. 

1.4 Literature Survey : 

A survey of some of the earlier work m the field of 
nonlinear discrete time systems will now be presented. As 
seen m the previous section, the nonlinearity in discrete 
time systems can occur either in polynomial form or in a 
piecewise linear form. This review is mainly concerned with 
the past work reported in nonlinear discrete time systems 
with polynomial type as well as saturation type of nonlinca- 
ritios . 
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Very little work was done on the analysis of discrete 
time systems prior to 1950. The early text books on Servo- 
mechanisms by Maccoll {.44], James et al [45] and Oldenbourg 
and Sortonus [46], all have some coverage on topics of 
sampled data systems, mostly limited to stability analysis 
and derivation of transfer functions. Since 1960 the field 
of discrete time systems has gone through a rapid change and 
the principle of optimal control, the concept of state varia- 
bles and stability studies by Lyapunov’s method [47, 48, 49] 
have been used for discrete time systems study. In stability 
studies the methods available for continuous time systems 
have been modified and applied to various types of nonlinear 
discrete time systems/ for instance Weaver and Leake [50] 
studied the stability and boundedness domain aspects of auto- 
nomous discrete time systems through Lyapunov's direct method. 
JDutchak and Sinitskiy [51] have formulated theorems to study 
the stability, asymptotic stability and instability of the 
solutions of nonautonomous nonlinear difference equations. 

Most of these theorems are very similar to Iyapunov* s theo- 
rems for continuous case. Liapunov's direct method has also 
been considered by Gordon [52] to study the stability prop- 
erties of perturbed difference equations and to analyse the 
stability aspects of discrete systems with periodic coeffi- 
cients in [53] . Two recent text books Freeman [54] and 
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CadzoW 1.55] written exclusively for discrete time systems 
discuss at length the linear system theory aspects of dis- 
crete time systems . 

As far as the analysis of difference equations are 
concerned, much less work has been done till recently. Milne 
and Thompson [31] and Jordan ct al [33] among others, have 
considered at length the various techniques that can be used 
to obtain the solution of linear difference equations with 
constant coefficients. A detailed description of various 
methods to analyse the linear discrete systems is also found 
in Freeman [54] and CadzoW [55] . 

A powerful mathematical tool devised for the analysis 
and design of discrete time system is the Z transform. 
Bagazzine and Zadeh [56] formally introduced Z-transform 
theory and applied it to automatic control systems. Barker 
[57] further extended this method and suggested what is now 
called the modified Z-transform method. The application of 
the modified Z-transform method to wide class of discrete 
time systems are discussed in the text books by Kuo [58] and 
Jury [59J. 

lindorff [60] has shown by introducing bilinear trans- 
formation that the Eouth Iiurwitz criterion and the Bode dia- 
gram can be applied to the linear discrete systems. It is 
also shown how the Byquist criterion and root locus technique 



19 


are valid techniques as applied to the discrete time model. 
Sedlar and Bekey [61] have applied newly formulated signal 
flow graphs to sampled data systems. Just as the methods 
available for solving a linear differential equation can be 
extended to linear difference equations, so also the methods 
developed for analysing nonlinear differential equations 
appear extendable to nonlinear difference equations. Details 
of applying the methods have to be altered somewhat, but 
much the same general approach appears possible. 

There is no general technique , however , to analyse the 
nonlinear difference equations as is also seen in the case 
of nonlinear differential equations. A considerable body of 
literature has been developed , however , on the subject of 
analysis of nonlinear differential equations. Analysis of 
nonlinear difference equations m the time domain has, how- 
ever, been studied to a much lesser degree than their counter- 
part, the nonlinear differential equations. 

Early work in this area can be traced to contributions 
by Trong [62] who presented a method for solving second 
order nonlinear difference equations containing small para- 
meter. The theory of first approximation which is widely 
used in solving nonlinear differential equations is extended 
to the field of nonlinear difference equations. Huston [63] 
has applied the discrete analog of the Krylov-Bogoliubov' s 
(K-B method) averaging method to the above equations. These 
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methods are very laborious and limited to autonomous systems. 
Transform techniques have been popular in control system con- 
text. Pai and Jury [64], Reddy and lagan [65], Iavi et al 
[66], among others, have applied respectively, convolution, 
multidimensional and modified multidimensional Z -transforms 
for the analysis of nonlinear difference equations. Non- 
linear time varying discrete systems have boon considered 
by Jagan and Desai [67] and analysed through multidimensional 
Z-transf orms. Here again only autonomous systems have been 
mentioned and these approaches are also fairly lengthy and 
involved. Alper [68] and Pu [69] have applied voltcrra 
series for the analysis of a class of nonlinear discrete 
systems and this approach is again limited to simple non- 
linearities. 

Aseltine and Nesbit [70] and Lindroff [71] proposed 
the incremental phase plane technique to analyse a non- 
linear sampled data system described by difference equations. 
Aseltine et al [70] have explained a method to obtain the 
difference equation to describe a sampled data system. A 
coordinate transformation m the incremental phaso plane is 
introduced [71] to identity int ersample oscillations which 
are undetectable in the original system. Much more about the 
methods is to be developed for example the possibility of 
their use to predict the nonlinear phenomena such as limit 
cycle behaviour. Pai [72] has proposed a new technique to 
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examine the necessary conditions for the existence of cer- 
tain limit cycles in nonlinear sampled data systems. This 
technique is based on the principle of specifying a certain 
repetitive output from the nonlinearity and finding conditions 
under which this sequence will be sustained. The analysis 
is based on a frequency dona in approach and requires an 
exhaustive search of all possible periodic modes. 

Second order digital filters described by nonlinear 
difference equations arc also considered in the present study. 
Such systems are analysed to investigate certain nonlinear 
phenomena such as limit cycle oscillations , jump phenomenon 
and super/sub-harmonic oscillations. We now survey some of 
the work done in the field of nonlinear digital filters. 

As mentioned earlier, because of finite wordlength 
available for the representation of signals, almost every 
digital filter is nonlinear. Effects of wordlength reduction 
in digital filters, such as limit cycles, overflow oscillations 
etc. have been known for a long time. Sinco tho early work 
of Gold and Ifcader [73], Kaiser [74], Jackson [75] and Ebert 
et al [76], m which these effects were mentioned for the 
first time, a large number of papers have appeared dealing 
with this subject. Methods have been reported for the ana- 
lysis of these effects, and measures for suppressing the 
unwanted phenomena have also been described, especially for 
second order wave digital filters [77, 78], The wave filters 
are related to classical networks and this relation transforms 
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the power function of the so called reference network into 
a function, known as pseudo power introduced by Fettweis 
[77]. The associated pseudo energy function can be used as 
a Lyapunov function for determining the stability of wave 
digital filters. An introduction to wave digital filters 
and their properties are also found m the above references. 

In the available literature on nonlinear phenomena due 
to wordlength reduction, a bulk of work has concentrated on 
limit cycles rather than on 3 amp phenomena and sustained non- 
linear oscillations. Basically there are two types of limit 
cycles, the first is due to overflow in tho adder when one’s 
or two's complement arithmetic is utilised and the second 
results from the rounding or truncation of the signal coeffi- 
cient product at the output of the multipliers. 

Four possible overflow nonlinear characteristics are 
reported in a recent paper by Clausen et al [79]. The study 
of overflow oscillations which can occur m the response of 
a secon " 1 order digital filter, whoso linear model is asymp- 
totically stable, was begun by Eberfc et al [76] and Sandberg 
[80], Ebert et al [76] have shown, for two's complement 
arithmetic, overflow oscillations can be sustained for appro- 
priately chosen initial conditions, while no input signal 
is present. These overflow oscillations are completely 
eliminated when a saturation kind of overflow nonlinearity 
is employed. However, if quantization effects are taken into 
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consideration, then the zero input response of a second 
order digital filter having an asymptotic stable linear 
nod el, can in fact possess overflow oscillations when 
employing saturation arithmetic. SairLberg [80] has shown 
that the amplitude of overflow oscillations can be made as 
small as possible by making the quantization effects snail, 
that is by using sufficiently largo number of bits in the 
representation of the data. 

Montgomery [81], Mitra et al [82-], Willson [83] and 
Claasen et al [79] among others studied the overflow osci- 
llations in second order digital filters with saturation 
arithmetic. Most of the work has been directed towards the 
question of existence of limit cycles, analysis [80, 84, 85], 
estimates on limit cycle amplitudes [86, 87, 88], stability 
criteria [89, 90, 91], and extension to higher order sections 
[92]. Parker and Hess [93] studied the limit cycle phenomena 
in a Somewhat greater breadth and developed new amplitude 
bounds based on Lyapunov' s method and a matrix oriented 
approach. 

Under forced situation, depending on the initial con- 
ditions two or more entirely different steady state responses 
can result for the same input signal [94]. A small change 
in amplitude or frequency of harmonic input signal may cause 
a jump between two distinct stationary solutions, usually 
known as jump phenomenon. Kristiansson [95] and Claasen 
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et al [96] have studied this effect m a forced second order 
digital filter with saturation overflow characteristics. As 
in continuous time systems, the response of a filter to a 
harmonic input signal, contains considerable subharmonic 
components. This phenomenon of subharmonic oscillation has 
been studied by Claasen et al [96] m a recent paper. Willson 
[97] has proposed a general error feedback circuit to minimise 
the presence of nonlinear phenomena, caused by the occurrence 
of adder overflow in the forced response of a second order 
filter. In general a large number of papers have already 
appeared in dealing with the subject of limit cycles in digital! 
filters and a comprehensive bibliography on the limit cycle 
problem is given by Kaiser [98] as well as Claasen et al [79]. 
On the other hand, there have been very few papers dealing 
with discussion of jump phenomena and super/sub-harmonic osci- 
llations in digital filters. 

1,5 Scope and Outline of the Thesis : 

The review of various methods to analyse the nonlinear 
difference equations arising m discrete time systems, given 
in the previous section has made it clear that there is no 
general way to approach a nonlinear difference equation. A 
few technique^ like Krylov-Bogoliubov' s averaging method are 
available to study the nonlinear discrete systems with cer- 
tain limitations. In this thesis a part of work is concerned 
with the application of a very recent scheme, namely multiple 
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time pcrturbational technique which has been extensively 
used for analysis of a class of differential as well as 
partial differential equations. This has been adapted to 
analyse nonlinear difference equations with polynomial type 
of nonlinearities. The main philosophy of this method is as 
follows. The solution x(k) is considered to be a function 
of two independent scales 1) a fast tunc rj and 2) a slow 
time t. The nonlinearity is essentially responsible for the 
modulation on a slow time scale of the solution obtained 
from the linear system (i.e. the nonlinearity assumed to 
be zero). The remaining work pertains to second order non- 
linear digital filters, which arc the basic building blocks 
for higher order filters, undar free and forced situations, 
limit cycle oscillations and nonet omcally decaying type of 
response, both outside the stability triangle under force 
free condition and the jump phenomena and subharmonic osci- 
llations under harmonically excited situation arc discussed. 

The chapterwiso summary of the work carried out in 
this thesis is given below : 

In Chapter two, a multiple scale pcrturbational tech- 
nique [99] has been developed and applied to linear as well 
as a class of nonlinear difference equations. For nonlinear 
equations the variation of output amplitude with input 
frequency namely the response characteristics are obtained 
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similar to those for continuous systems [5]. A nonlinear 
difference equation of Van dejr .Pol type has been analysed 
and qualitative inf orraation such as limit cycle oscilla- 
tions is obtained. The stability of such limit cycle 
oscillations is studied through a vanatiunal technique. 

The scheme [99] used for the above analysis is obtained 
through the contral difference operator property. This 
scheme is an improved version of the one reported by 
Lrishnan and Subranaman [100J, using the forward difference 
operator property, to analyse a class of difference equa- 
tions. Both arc suitable for qualitative kind of analysis 
such as limit cycles, oump phenomena etc. under steady state 
condition, whereas the scheme reported m [100] is not 
satisfactory for quantitative analysis to obtain the transi- 
ent solution. Whenever a transient solution of a given 
system is required, a modified definition to the above 
schemes is given, which can be applied to obtain the appro- 
ximate solution closer to the exact solution. The modified 
scheme is formulated out of experience and the mathematics 
involved is less than the proposed method in [99]. In 
practice most of the nonlinear systems require only quali- 
tative analysis due to complexity m the nonlinear part of 
the system. Examples are given to illustrate tho modified 
technique and tho results arc compared with those obtained by 
another recent multiple scale perturbatiaval technique [101] 
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which is based on Taylor series expansion. 

Chapter three deals with the important concept of 
super/sub-harnonic oscillations in nonlinear discrete time 
systems under strongly excited situation. The discrete mul- 
tiple time perturbational technique which has been discussed 
m detail in Chapter two as well as the well known technique 
of 'Harmonic Balance' method are adapted m the present 
study. A Duffing kind of equation is considered and the 
possible super, sub and ultra sub-harmonic oscillations arc 
investigated. The results deduced from the proposed methods 
are compared with the ones obtained by computer simulation 
of given system equation. 

Chapter four concerns with study of second order digi- 
tal filters with overflow saturation nonlinearity under 
zero input situation. The survey of the past work in this 
area reveals that there arc no limit cycle oscillations with 
this typo of overflow characteristics, when the quantization 
effects arc ignored. However, zero input limit cycles do 
exist outside the stability triangle with saturation arith- 
metic. Three methods have boon proposed to delineate the 
regions of different limit cycles in the a-b parameter plane 
outside the stability triangle. Another important and 
interesting contribution hero is the identification of a 
region outside the stability triangle, in which the linear 
system is unstable and the nonlinear system is asymptotically 
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stable for a particular set of initial conditions. Tho 
region in the a-b plane and the set of values m the initial 
condition plane for which the aforesaid behaviour is possible 
are indicated. The system equation is simulated and the 
results are compared with the theoretical findings. A 
brief description of these two aspects namely the limit 
cycle oscillations and the monotonically decaying nature of 
the response, outside the stability triangle* m the a-b 
parameter plane is explained in [102]. 

Chapter five considers forced second order nonlinoar 
digital filters. The contributions m the field of jump 
phenomenon are limited and only two papers by Kristiansson 
[95] and Claasen et al [96] are available in literature. 

No mention is made anywhere in these references, of the 
range of the input amplitude for which there exists a jump 
behaviour or nonlinear sustained oscillations for a speci- 
fied set of filter parameters [a,b]. A new method is propo- 
sed to find the range of input amplitude for which there 
exists always a jump [105 ] t for a specified frequency and 
the filter coefficients. This is a simple graphical constru- 
ction of derived conditions for a specified frequency and it 
gives regions m the a-b plane in which either jump or 
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nonlinear sustained oscillations arc possible. The discu- 
ssions on subharmonic oscillations, methods adopted to 
eliminate the unwanted nonlinear phenomena and the verifica- 
tion of results by simulate jn form the renaming part of 
this chapter. 



CHAPTER 2 


FREE AND WEAKLY FORCED NONLINEAR DIFFERENCE 

EQUATIONS 


2.1 Introduction ; 

The present chapter is concerned with nonlinear difference 
equations with a polynomial kind of nonlinearity. With the 
advent of high speed digital computers, discrete time systems 
described by difference equations have been an important 
and practical model for the scientists and engineers, especia- 
lly the control engineers. As a result, the theory of 
difference equations has been rapidly developing in the 
literature. Much of this research has been devoted to the 
development of results which parallel those for differential 
equations. 

In the theory of oscillations, the method of linearisa- 
tion has acquired a special significance and is widely used 
to analyse nonlinear oscillatory systems. This subject is 
covered in many text books including krylov-Bogoliubov J_12] 
who formulated the powerful scheme presently known as Krylov- 
Bogoliubov (K-B method) averaging technique to analyse a 
large class of nonlinear continuous systems. An amplified 
version of this method known as Krylov-Bogoliubov - Mitropolski 
(KBM) method dealing with the asymptotic theory of oscillatory 
behaviour of a class of nonlinear systems is due to Bogoliubov 
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and Mitropolski [13]. The K-B averaging method has also been 
adopted by Huston [65] and Trong [62] to analyse a class of 
nonlinear discrete time systems. The application of K-B 
averaging technique to nonlinear difference equations has 
however been limited to autonomous a ses and the mathematical 
computation involved is lengthy. Transform methods [64-67] 
nave been used extensively to solve linear and nonlinear 
difference equations. These methods, again are restricted 
to systems without forcing. 

A recent method, namely the multiple time perturbation 
technique, introduced by Cole and Kevorkian [14], has been 
applied extensively to study a class of second order differen- 
tial and partial differential equations by a number of workers 
including Cole [16] and Nayfeh [18]. The same scheme has 
been extended by Tiwari and Subramanian [21] to study the 
response characteristics in weakly third order nonlinear 
differential equations under the free and forced situations. 
As far as analysis of nonlinear difference equations are 
concerned, this multiple time scaling approach does not seem 
to have been attempted so far except for a recent work by 
rloppensteadt and Mira inker [101] who considered a general nth 
order discrete time system described in matrix form. In 
[101] the discrete time model has been obtained through a 
Taylor series expansion. A new technique is adopted in this 
thesis to develop the two discrete time scaling model through 
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the known properties of finite difference operators. A 
brief review of the properties of the finite difference 
operators is given in the next section. The proposed two 
time scaling approach is quite general and can be applied to 
study linear as well as a class of lonlinear difference 
equations with and without a forcing function. To illustrate 
the idea, a Van dor Pol type of problem is considered and 
limit cycle oscillations under steady stabe concitions are 
obtained. The stability aspects of such limit cycle oscilla- 
tions are also analysed through a variational approach. A 
discrete version of Duffing equation is ti en studied fully 
under both free and weakly forc-.d solutions. The discrete 
model of continuous time Duffing equation is ob. lined by 
exploiting some properties of t ie central difference opera- 
tor. The variation of steady state amplitude with input 
frequency for weakly forced situation, known as response 
curves are derived and analysed. 

The proposed discrete two time expansion method using 
finite difference operator properties is more suitable for 
qualitative analysis of nonlinear phenomena, namely descri- 
ption of the steady state solution like limit cycle oscilla- 
tions, response characteristic curves etc. On the other hand, 
the transient solution results deviate somewhat from the 
exact solution. To obtain more accurate result closer to 
the exact solution, the initially proposed scheme of two 
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time scaling is redefined and the results obtained using the 
second form are compared with those obtained using the method 
proposed in [101], as well as the exact solution obtained bjr 
computet- simulation of the £iven system equation. The modi- 
fied scheme proves to be a better a pproach for discrete 
problem analysis under certain conditions which are explained 
later on in the analysis. The results obtained are verified 
by computer simulation. 

2.2 Multiple Time Perturbation Method in Continuous Time Systc 

Before proceeding with the basic analysis it is useful to 
review the method of the multiple variable expansion procedure 
employed to analyse a class of second order nonlinear differen- 
tial equations. The two time scaling method , a special case 
of multiple scale perturbational technique has been used to 
analyse a variety of weakly nonlinear differential equations 
(sometimes known as small perturbation problems). The two 
time scaling technique assumes variation of the dependent 
variable x interms of two independent time scales, a fast 
time and a slow time x, The motivation for introducing the 
two independent variables has been explained m the introdu- 
ctory chapter. In general both these independent time scales 
are assumed to be a function of a small parameter p., usually 
known as perturbation parameter and are expressed as [20,104] 
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2 3 

T) = t(l +w^i -Mo + . . ..) the fast tune scale. 

, 2 3 ( 2 . 1 ) 

t - t(p + S^i + S^i + . ..») the slow tune scale 

where the constants w and S are to be determined so that 
the asymptotic expansion is uniformly valid. The constants 
are included in the definition of slow time scale to allow 
for general dependence of nonlinear function on \i [104]. 

Cole [16] and Uayfeh [18] have considered the definition 
of the slow variable in a slightly different manner by assu- 
ming the constants = 0. However # as far as the nonlinear 
systems analysis is concerned, second approximation is possi- 
ble in theory and impossible m certain situations. Under 
such situations the definition of independent scales can be 
simplified by assuming ^ = a n = 0, that is 



Generally the application of multiple scale procedure 

is limited to weakly nonl- near systems and the contributions 
2 3 

of terms }i , y are less significant and hence for the first 
approximation the independent scales defined in eqn. (2.2) 
can be used for analysis. 

Despite the fact that solutions generated by the two 
time scaling method appear to satisfy the equations uniformly 
for all time, in general it is known that the approximate 
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solution is a valid approximation only for times 0(l/p) as 
|i approaches zero. The relation among the dependent variable 
x, the perturbation parameter \i and the two time scales is 
given as in asymptotic serie" as follows. 

2 

x(t) = x(r),x) = x 0 (t},T) + px^(r), v) + \i x 2 (n,v) + .. (2.3) 

where \i is small, that is 0 < p « 1.0. In eqn. (2.3) 

x 0 (ri,x) is known as 'generating' or 'base' solution and 
x-j_(r),x), x 2 (r),x), etc. are the first, second, etc. correction 
terms. It is to be noted that for weakly nonlinear syst sms , 
the terms upto first correction term in the series solution 
gives fairly accurate approximation for all time. The accu- 
racy of the solution can be further improved with the addi- 
tion of more correction terms. For practical applications 
terms upto x^(r),x) are sufficient. 

The basic idea behind the two time expansion technique 
may be illustrated by considering a general second order 
differential equation of the term 

x + |i f (x,x) = 0 (2.4) 

where a> 0 is the natural frequency, \i the perturbation 
parameter and f is a nonlinear function of indicated variables. 
It is assumed that the function f is differentiable with res- 
pect to its arguments. ' ♦' denotes the derivative with res- 
peot to time t. 
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The solution x(t) to the egn. (2.4) is considered as 
a function of two time scales and is given by 

x(t) = x(t),t) (2.5) 

which in turn can be expressed in an asymptotic series form 
as 

x(ti,t) = x 0 (t),t) + |ix 1 (ti,t) + |i 2 x 2 (t),t) + .. (2.6) 


Now the dependent variable x is a function of two 
independent variables n and t, the time derivatives x and x 
can be expressed as partial derivatives with respect to die 
two time scales as follows : 


* = = asinxil = (1 + «, i(1 2 + w 2 |j5 + 

2 ^ 

4- (\i 4- S^Jl ^ + ♦ ••) 
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Then the given eqn. (2.4) becomes 


+ 2}i + I» 2 L2-1 - 3 - Z? V - 

3 n 3 1 d 3 r) 


+ 2s, l -Si lLi -T i + 

1 311 3 t 


x(n ,t 


h { f( V 57r } 


+ ' llXl x-x + 


( afi + li 

3t) 3 t 


■)| 1 } 

x=x„ 


combining the eqn, (2.6) with the above expression. and 

o 2L 2 

collecting coefficients of p , \i , p , ... terms and equa- 
ting them to zero separately, the following set of linear 
partial differential equations are obtained 


a z 0 (n,v) 2 

2 +w o X 0^ ,T ^ ■■ 0 

2 *** 

3 X 1 . 2 _ ^ 3 2 x^(ti.t) 3X 0n 

— T- + “ O = - 2 • g - Vv • - f (x 0» rn ) 

3 T) 3 TQ 3 V 


( 2 . 8 ) 


(2.9) 


9 X 2 (t),t) 2 


2 — + w o 


3 ^(r),*) 
3 1 }t 


- [2 -! 3 2 * ° (T) ’ T) + 2s 3 Vl'ji + 
rts 1 3 n 9 1 


l*lf§ 


a x, a x„ "3 f 

+ ( 3~ + a“ } ^ x I - 

J 3 x~ 

O 0 

* 3_fo • 3 x 0 

a t) x ~ a ri 


(2.10) 
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Th.e solution to (2.8) is the base solution and is expressed 

2t S 

= A(t) cos a) o n + I(t) sinw Q n (2.11) 

where A(x) and B(t) are the slowly i varying amplitude functions 
of slow time variable t and can be evaluated by considering 
the first correction equation given m (2.9). 

To obtain the first correction terra, the generating 
solution is now substituted m the right iiand side of the eqn 
(2.9) and supression of secular terms by collecting coeffi- 
cients of cosw q t) and sin w q t} terms and equating them 
separately to zero result in coupled first order differential 
equation in A(x) and B(x). The solution to these coupled 
equation gives the explicit expressions for the amplitude 
functions A(t) and B(x). It is to be noted here that the 
explicit solutions for the amplitude functions A(x) and 
B(x) are possible for linear and for some simple nonlinear 
cases. However^ numerical methods can bo adopted to extract 
the information regarding tue amplitude functions. Then the 
eqn. (2.9) can be solved for x^(r),x) with nonsecular terms 
on the right hand side. Then the approximate solution upto 
first correction term, that is upto order n, is given by 

x(k) = x(t],t) * x 0 (p,t) + (i acj/*),*) . (2.12) 
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In. a similar way the second correction term, lamely 

03,11 a<i< i e ^ with the first corrected soli bion given 

m eqn. (2.12) to improve the accur .cy upto order \i . This 
kind of repeated addition c ~ correc >ion terms are possible 
in theory and generally impossible m practice r 

An excellent account of this method of stily, both in 
depth and breadth, can be found m Uole [16], lu.gteh [18 J 
and Reiss [105]. 

Inspite of many advantages, the multiple tir e expaa, .ion 
procedure suffers due to nonavailability of general procedure 
to express the two independent time scales in. terns of tie 
perturbation parameter p. However, in a very recent paper 
bovine and Obi [106] have analysed nonlinear differential 
equations adopting the two time scales generated from the 
given differential equation. 

Preliminaries in the finite difference schemes which 
are useful to derive the multiple discrete time perturbation 
technique are now given. 

2.3 Finite Difference Schemes ; 

Before proceeding with the derivation of discrete multi- 
ple time perturbational technique, it is useful to review some 
standard results in finite dif-J orenco schc >os. These results 
will be used further on m the analysis, for a given real 
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function of real variable f(x), following first and second 
difference operators are defined |_33J. 

(i) forward difference operator ( A ) 

A f(x) = |.f (x+h) - f(,x )J 
A^f(x) = [f(x+2h.) - 2f(x+h) + " (x)] 

the step size or increment h can be normalised t,o unity withbut 
loss of generality and the above equations take the form 

A f ( x ) = f(x+l) - f(x) j 

A 2 f (x) - f (x+2) - 2f (x+1) * f (x) (2.13) 

(ii) central difference operator (6) 

6f(x) = f(x+-£-) - f(x^-J-) 

l 

6 2 f(x)= f ( x+1 ) - 2f(x) + f(x-l) (2.14) 

the first central difference operator 6 may also be defined 
as follows « 

6f(x)= [f (x+1) - f(x-l)]/2 (2.15) 

(ni) backward difference operator (ST ) 

<?f (x) = f(x) - f ( x+1 ) 

V^f(x)= f(x) - 2f(x+l) + f(x+2) (2.16) 

Then generalising to a function of twe variables x and 
y the partial difference operators are defined as follows : 
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(iv) Forward partial difference operators j 

A = f(x+l ,y) - f (x,y) 

A y f(x,y) - f (x,y+l) - f(x,y) 

A x f(x,y) = f(x+2,y) - 2f(x+l,y) + f(x,y) 

A yf(x,y) = f (x,y+2) - 2f (x,y+l) + f(x,y) 

A x Ay f ^ x,y ) = f ( x+1 * y +1 ) - f(x+l,y) - f(x,y+l) + f(x,y) 

, \ (2.17) 

(v) the central partial difference operators : 

6 x f(x,y) = f(x+^,y) - f(x-4,y) 

6yf(x,y) = f (x,y+4r) - f(x,y-£-) 

6^f(x,y) = f (x+l,y ) - 2f(x,y) + f(x-l,y) 

6yf (x,y ) = f(x,y+l) - 2£(x,y) + f(x,y-l) 

«5 x 6 y f(x,y) = f(x+4-,y-40 - f(x+i,y-^) - f(x-4 » y+i) + 

±(x-%, y-£0 „ (2.18) 

In the following section, the development of discrete two 
time scaling technique using the above finite difference 
operator properties will be illustrated. 

t 

2.4 Discrete two Tune Scales Expansion : 

Knowing the properties of various finite difference 
operators discussed in the previous section, a two time 
scaling model can easily be obtained for discrete time 
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systems. The aim is to derive a stable (stable in the sense 
that the solution obtained by the derived technique is same 
as that of exact solution qualitatively) discrete multiple 
time perourbation technique, which when applied to difference 
equations would give approximate solation for all discrete 
time instants. Roscue [107] formulated a new stable difference 
scheme for differential equations which have exact solutions 
m certain simple situations. Numerical solution of partial 
differential equations has been studied by Smith [108] by 
replacing the derivatives with valid difference schemes. 

Based on these schemes a valid discrete time scaling scheme, 
analog to two variable expansion method in continuous time 
systems, is derived as follows. 

let T be a general difference operator representing 
a forward, a backward or a central difference operator. 

Let I denote the set of nonnegative integers. The two 
independent time variables are defined as follows : 

2 *3 

r) = lc(l u 2 (1 *•*) fast time scale 

0 *7 

x = k(l + s-jji + s 2 iv + slow time scale (2.19) 

where kel, w and arc constants to be determined so that 
the asymptotic expansion is valid uniformly. Note that the 
general dependence of nonlinear function upon \i demanding the 
inclusion of the constants m the definition of x [104] is 
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not true in the discrete time definition. This will be 
illustrated later on m the analysis. 

Now the given dependent variable x(k) is a function of 
two independent scales r) and x and is expressed as 

x(k) = x(p,x) (2.20) 

As mentioned earlier T is a general difference operator, 
which operates on the associated function with respect to 
k, operating the eqn. (2.20) with T 

Tx(k) = Tx(p,x), 

In the above equation, the arguments of the variable x are 
p , x which are m turn functions of k. Under this situation 
the right hand side operation of the above equation is 
defined as follows 

Tx(k) = T x(p,x) Tp + I x(p,x) Tx (2.21) 

T) T 

whore and T^ are the general partial difference operators. 

It is to be noted that irrespective of the kind of opera- 

2 3 

tion that T exeoutes, Tp = (1 + w -jji + co + ....) and 

f 2 "5 

Tx = (|i + s^|i +■ + .....). This invariant property can 

easily be checked as follows for various finite difference 
operators. 

(i) T = A , the forward difference operator. 

2 ^ 

Tp = A -rj = A (1 + 2 ^ + . . . ) Jr 

2 *2 2 3 

= (1 + + to 2s + * * • ♦ ) ^ ^ = (1 + + a> 2 ^ • • * ) 
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2 3 

Tx = A t = A (y + s 1 (i + s 2 \r + . . .)k 

2 3 2 3 

= ( \i + s + s 2 ^ + . • • ) A k = ( p + s-jjx + s 2 (x t • • • ) 

(in.) T =v t the backward difference operator 
Hr) = v h = (1 + oj-l^ 2 + + . . . ) k 

= (i + W j\l^ + w + •••#) 7 i 

= (1 + W-jji 2 + w 2 ,i 3 + ....) 

2 3 

Tt = sjx = <7([i + s- L n + s 2 H + . . . ) k 

2 3 2 3 

= ( \i + s j\ l + s + • . • ) ^7 lc = ([i+ s lJ l ^ s 2^ x ^ • * * , 

(in) T = 6, the central difference operator 

Tr| = 6r] = 6(1 +w -j_^ 2 + a) 2 ^ + . . . ) k 
2 3 

= (1 + w ^ + w 2 ^ + * • • ) 6 fc 

= ( 1 +■ t w 2 ^ « • • ) ( k + -£■ — k + 4r ) 

2 3 

— (l + ^ -j^|i + ,- u. 2 (i + • »•) 

2 5 2 5 

Tx = 6t = 6(p + s 2 p + •••)k = ([i + s^x + s 2 }r + ..) 

2 5 

= (n + s-jjx + s 2 P + • ...) 


(iv) 1 = 6, the central difference operator defined as m 
eqn. (2.15). 
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O 7L 

Tr) = 6r) = 6(1 + w + w 2 p^ + , . . . ) 

= (1 + u ^|i 2 + <»> 2 ^ + . ..)6k 

= (1 + to j\i + w 2 ^ + • * * ) 

= (1 +« lt i 2 + « 2 p 5 + . ..) 


2 3 

Tt = 6 % - 6(p + s-^p + s 2 p. + . ,.)k 

2 *5 

= (n + s^n + s 2 n + . ...)6k 
2 *5 

= ( (1 + S^Jl + S 2 (i + ....). 


Based on the above invariant property eqn. (2.21) can be 
rewritten as follows 


Tx(k) = (1 + «-,p 2 + a) 0 \i? + ...) T x(r ) ,%) + 

JL t - YJ 

+ (p +■ s^p + s 2 p^ + ....) T x(t)jt) . (2.22) 

Operating on the above equation again by T, wo get 

f 2 x( k) = (1 + “-jji 2 + w 2 p 5 + . . . ) 2 T 2 x(p,v) 

+ 2(1 + w.^p 2 + w 2 p^ + « . . . ) x 

(p + s^i 2 + s 2 p 3 + ...) H T % x(t),t) 

+ (p + s-^p 2 + s.p^ + ...) 2 T 2 x(t],v) , (2.23) 

Eqns. (2.22) and (2.23) can conveniently be written as follows 


Tx(k) = I x(r),v) + pT x(ti,t) + 

+ p 2 [ T x(t),t) + s 1 T x x(t),x)] + 0(p 5 ) . (2.24) 
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T 2 x(k) = l 2 + 2^1^. 

+ H 2 [2 x(t},t) + 2s 1 T ri T T x(r],x) + T 2 x(ti,t)] + 

0(\i 3 ) . (2.25) 

Equations (2.24) and (2.25) give the general discrete two 
time expansion scheme. A stable difference scheme can now be 
obtained replacing T by & in eqn. (2.24) and T by 6 in eqn. 
(2.25). It is to be noted that the replacement of 1 by 6 
m both the equations gives an improved stable scheme, but 
this scheme is not accessible for analysis of the resulting 
difference equation, both by analytical means as well as by 
computer due to the noninteger arguments appearing as sub- 
scripted variables. Thus, the scheme used in this thesis for 
the analysis of discrete time system is given below : 

Ax(k) = A x(t),t) + (i a t x(t},t) + 4 L o>iA ^ x(tj,t) 

+ A T x(ti,t)] + 0(p 5 ) 

6 2 x(k) = 6 2 x(rj , x ) + 2\i 6^ x(p,t) + p 2 L 2a) i^ 2 x(n,T) 

+ 2s 1 6^ 6^ x(t),t) + 6 2 x(r},v)] + O^i 5 ) (2.26) 

where 

A is the forward difference operator, defined as 


A x(k) = x(k+l) - x(k) and 
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6 is the central difference operator, defined as 
6 x(k) = x(k+^-) - x(k— £■) 

and a ^ a x > 6^ and 6^ are the partial difference operators. 

In general, most of the discrete time systems are des- 
cribed by difference equations of the following form 

x(k+l) + ax(k) + px(k-l) + g[x(k), x (k-1)] = 0 
rather than using difference operators, for instance 

6 2 x(k) + a6x(k) + px(k) + g[x(k), 6x(k)] = 0 

and hence the expressions given in eqn. (2.26) are expanded 
and rewritten as follows ; 

x(k+l) = x(t}+1,t) + p|x(r) ,t+1) - x(t),t)] 

+ (i L u> ,x(tj+1,t;) + s-^x(r) ,t+1) 

- ( « + s i^ + O(n^) (2.27) 

x(k+l) + x(k-l) = x(t]+1,t) + x(ti-1,t) + 

- 2x(ti+^,T— i-) - 2x(r)~^,TT+^) + 2x(r}-4r,'C-£)] 

+ \ i 2 { 2 Ml z(n+1,T) - (4^+2) x(p,t) + 2w 1 x(t}— 1,t) 

+■ x(t|,t+1)+x(ti,t~ 1) + 2s 1 (_x('n+4,T+4 r ) 

- x(t]+-£,t~£-) - x(t}-£,t+£) + x(rHr,T“ 4 )]} + 0 (^ 5 ) (2.28) 
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from the above two equations 

k(k-l) = x(t]-1,t) + |i[_ 2 x(t)+-3t, t+-^-) - - 2x(r)~ £-,?+•£•) + 

2x(n-£,T-£-) - x(t),t+ 1) + x(r),x)] + 

[ ? { U) -jE^tl,*) “ (5 U1 “ S.j + 2) + 

2 a) ^x(t)-1,t) + (l~s^) x(r),T+l) + x(r),T-l) + 
2s 1 [x(n+^-,T+4-) - x(r)+-^,T-4r) - x(r)-4, T+-J-) + 

x(t)-4-,t~-^) ] } + O(^) . (2.29) 

The equations (2.27) - (2.29) are used to transform the 
given difference equation into two discrete time model. The 
system equation and its analysis using the above equations are 
described in the following sections. 

Remark : 

The operator T m eqn. (2.24) can also be replaced by 6, 
defined in eqn. (2.15) as 

6 x(lc ) = (*) 

but the difficulty here is the order of the system is increased 
by one. So, whenever 6 is appearing by itself in any equation, 
the replacement is by the forward difference operator A , 
namely 

6x(k) = A x(k) = x(k+l) - x(k) . 
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On the other hand when 6 is appearing along with higher 

2 3 

order difference operators nanely 6 ,6 .... etc., the 
replacement for 6 is as defined in (*). However , the replace- 
ment of 6 with usual central difference operation namely 

<5 x(k) = x(k+£) - x(k—i !•) 

is possible when explicit expression for x(k) is known. 

2.5 System Model ; 

A class of nonlinear discrete time systems considered 
in this chapter for analysis is a second order difference 
equation of the form 

x(k+l) + ax(k) + px(k-l) + pyg[x(k), x(k-l)] = \i F(k) 

( 2 . 30 ) 

where kel, is the indexing discrete time parameter, a,p,y 
are constants and \i is a small positive constant. F(k) 
is the external input, usually a periodic function. 
g(*, •) is a general nonlinear function with indicated varia- 
bles . 

For = 0, eqn. (2.30) reduces to a linear difference 
equation and is given by 

x(k+l) + ax(k) + fJx(k-l) = 0 „ (2.31) 

It is assumed that a,{3 are so chosen so as to guarantee 
a bounded response. The different type of solutions for this 
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linear case are given in Appendix A. It is to be noted 
that for [a,{3] lying anywhere within the stability triangle 
ACG shown m Fig. 2.1, a bounded response (that is decays 
to zero f olution) is obtaini led. For [a,p] lying outside 
the stability triangle including the points A and C a *e the 
triangle aCG, an unbounded response is obtained. On the other 
hand an oscillatory solution is possible on the line AC and 
at point G (except at points A and C), and steady state 
oscillatory solutions are possible on the boundary lines AG 
and CG. With this knowledge of the properties of the stabi- 
lity triangle, the analysis of the nonlinear system (for 
p. 4 0) will be provided in the next section. 

2.6 Analysis : 

With the complete knowledge about the solution of the 
linear difference equation (2.31), the approximate solution 
of nonlinear equation (2.30) for \i 4 0, will now be obtained. 

Before proceeding with the detailed analysis, it is con- 
venient tro make a few commer fcs regarding the nature of the 
solution for the linear system (for \i = 0). An observation 
of the types of bounded response reveals that for [a, (3] 
lying on the boundary of the stability triangle ACG (points 
A,C excluded), a steady state periodic solution is obtained. 

For lying m the inferior of ACG it is convenient to 

define a detuning parameter N such that for the force tree 

I 1 * * 

586 ^| 
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system the discript ion takes the term 

x(k+l) + ax(k) + x(k-l) + pNx(k-l) + pyg[x(k) ,x(k-l)] = 0 

( 2 . 32 ) 

where IT is the detuning parameter defined as 

« (P - 1) . 

The advantage gained by such a transposition of terms 
is the generation of a ’base' solution as a periodic one 
(corresponding p° term in (2.32)), For detuning parameter IT 
of the order of unity or less, a good approximation to the 
true solution is then easily obtained as pointed out later in 
the analysis. Again for the forced situation, that is for 
F(k) ^ 0, it is likewise convenient to define smother detuning 
factor that takes into account the deviation between input 
frequency #uo(P(k) assumed periodic) and linear system fre- 
quency. It can be expected that for small values of this 
detuning parameter a strong (large amplitude) response will 
then be obtained, 

Considering F(k) = Pcos w k and adding the term 
2 coso) x(k) on both sides of eqn. (2.30) and rearranging, we 
get 

x(k+l) - 2 cos w x(k) + x(k-l) + plTx(k-l) + pyg [x(k),x(k~l)] + 


pM x(k) = pF cos ca k 


(2.33) 
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from e^ns. (2.34-) and (2.35), the following sets of linear 
difference equations are now obtained by considering the terms 

_ -j o 

with coefficients [i , (i , \i f ...... and equating them indivi- 
dually t ' zero. 

|i° terms ; 

x 0 (r)+l,t) - 2 cos w x o (ti,t) + x 0 (t)-1,t) - 0 (2.36) 

\l~ terms ; 

x 1 (ti+1,t) - 2 cosa) x- L (r],'c) + x 1 (ri-3.,i:) = 

- { 2x 0 (TJ+-J-, T+£) - 2 X Q (riH»H) - 

2x q (ti-4-,t+^) + 2x 0 (r)-4-,T-i-) + Nx 0 (t)-1,t) 4 
YS|> 0 (ii,t), x 0 (ti-1,t)] 4- Mx 0 (ti,t)} 4- 

3? eosw r) (2.37) 

|i ^ terms : 

x 2 (ri+l,'c)-2cos a ) :£ 2 (t),t)4-x 2 (ti-1,x) = - 2x 1 (t)+£,t4-|-) - 

2x 1 (ti+^,'T-4-) + 

2x^(rj — Jr»x4-4r) 4- 2 x-^(t) \ , x—4r ) 4- Nx-^(t)-1,t) 4- 

p 

+ coefficient of n in g 4- 2“-^ 
x 0 (r}4-l,i:) - (4 U 1 4- 2) x q (ti,'c) 4- 2^ 1 x 0 (t)-1 i t) + 

x 0 (n,T4-l) 4- x 0 (ti,t- 1) + 2(» 1 4N)[x 0 (n4-J-,x4-4) - 

x 0 (r)+4-,x-4) - + r 0 (TT4» T -4)] } 


(2.38) 
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The solution to eqn. ( 2 . 36 ) in the generating or base 
solution and is obtained as follows : 

The characteristic equation is 

2 

m - 2 cos w ei + 1 = 0 
then the roots are 

2 2 
2 cos a) + f4 cos w -• 4 2 cosu + 3 yA -4 cos w 

m ± t 2 = 5 2 

= cose*) + 3 sm w 
Magnitude R 

Phase 0 = tan "" 1 = 

therefore the solution is 

XqCt'Ij'c) = A(t) cosu r) + B(t) sinu rj . (2.39) 

A(t) and B(x) are constants for linear case, that is 
for p = 0 and are otherwise slowly varying amplitude function 
depending on the slow time variable x. It is of interest to 
obtain explicit expressions of the amplitude variations when 
the system is weakly nonlinear, that is |i ^ 0. The method 
of obtaining the amplitude variations A(x) and B(v) by su- 
pressing secular terms generated on tho right hand side of 
eqn. ( 2 . 37 ) proceeds as follows : 

Lot the difference interval be ad 3 usted such that there 
exists some integer L such that 


= V cos^»> + sin^ w = 1.0 
— tan "" 1 (tan co ) = w 
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0 ) - 


271 


Expanding the nonlinear function g(*,*) in eqn. (2.37) in a 
discrete fourier series [109] as : 


P I / 2 

o . I n 2ukn 


SLx 0 (t),t), x 0 (t)-1,t)] = + * C n cos + 


where 


1/2 

l \ Sln 

n=l 


2-rckn 

1 


(2.40) 


"n 


It 

I 

g( *, 0 



T)=l 




L 



2 

“ I 

l 

g( *2 *) 

ooa Ssia 

T)=l 




1 



OJfcq 

1! 

l 

§(*>•) 

sin , n = 1,2,3**. .L/ 2 

T)=l 



and 

b-, are the 

coefficients of cos y 


sin respectively and are given by 

1 


3 1 “ I ^ «(*»•) coS 

T)=l 

1 
r 

g ( • , * ) sin ^Sp 4 (2.41) 

T|=l 

Then eqn. (2.40) can he expressed with c Q , c-^ and b^ as 

c. 


h -s 


2 l 


gl_x o (T),T), x 0 (ti-1,t)J = ^ ^ cos + b ± sin 




= v>r + c-^ cos w n + b-^ sinoj t] „ (2.42) 
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Prom eqn. (2.59) 

0) (i3 

x (t^, x+%) s [a(t+£) cos 2 + B(t+4-) Sin 5] cos w rj 

- [A(vH-) sm ^ - B(t-kJ-) cosw /2] sinmr] 

X 0 ( T )H> ''Hr) - L A ( X H) cos ^ + B(t- 4) sinSj] cos w q 

- [A(t—|-) smm/2 - B(t— £-) cos u /2] sinwr] 

x 0 (ri-4-, T+ i) = L A ( T+ D coS m/2 ~ b ( t H) s:L nm / 2] cosm q 

+ [A(th-£) sin m/2 + B(t+£) cos m /2] sin m n 

x 0 (bH» ''*-£•) = L a ( t H") cos m /2 - B(t—£) sinm /2] cosm q 
+• [A(t— £-) sin 05 /2 + B(t— £■ cos u / 2] sin w t) 

x 0 (q-l,x) = [A(t) cosm - B(x) sinm ] cosm q 

+ |_A(t) sm m + B(x) cosm ] sm m q < (2.45) 

Substituting title eqns. (2.42) and (2.45) in (2.57) and 
supressing the secular terns by equating the coefficients of 
cos m q and sinmq to zero separately, the following equations 
for amplitude variations A(t) and B(x) are obtained : 

6 A(x) = 4-5—0572 LMBC*') - \ + K(A(x) sinw + B(t) cosm )] 

6 B(x) = X~S ln^ "]2 E F + °i ~ m ( t ) ~ N(A(t) cosm - B(x) 

smm )] # (2.44) 
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Explicit expressions for A(t) and B(t) can be obtained 
from the above coupled equations if the system is linear or 
if the nonlinearity is of simple form. Such explicit expre- 
ssions ive however, not generally possible in most of the 
cases, numerical techniques may then be resorted to get 
the required result. The solution to eqn. (2.44) and sub- 
sequent substitution in eqn, (2.39) provides the first order 
approximate solution and is given by 

x(k) * x o (t),i;) = A(t) cosw t) + B(t) sinw rj , 

Then the eqn. (2.37) is solved with the remaining non- 
secular terms on the right hand side of the equation. Let 
the solution be 

Xi( t j,t) = C(t) cosw r) + D(t) sinw p + steady state solution 

due to noneecular terms 

The amplitude variations a(r) and D(t) can be determined 
considering the eqn. (2.38) and supressing the secular terms 
arising in the right hand side due to substitution of x g (ri»' , r) 
and x^tutO* Knowing the complete solution of now 

the approximate solution is given by 

x(k) * x q (t),t) + |i x 1 (ri, t) . 

A similar procedure can bo followed to include the 
higher order correction terms namely 3 C^(ti#t) etc. 
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But tlie addition of the correction terms with the generating 
solution is possible when the system is linear or when the 
system is described with simple nonlineanties . The inclu- 
sion of first two correction- terms, that is the approximate 

2 

solution to the given equation upto the order of ji is illu- 
strated by considering a linear oscillator in example 2.7.1. 

The result obtained shows that the improved approximations 
are of course possible by considering higher order terms m 
the assumed series form. 

The next section considers various examples, linear as 
well as nonlinear cases which will illustrate in detail the 
analysis technique as ewell as the applicability of the two 
variable expansion scheme to a clss of difference equations. 

2.7 numerical Examples : 

The suitability of the proposed multiple discrete time 
perturbation technique to a general class of difference equa- 
tions is illustrated in this section through some numerical 
examples, linear and nonli tear equations are analysed and 
the results are compared with the exact solution of the given 
equation (simulated as a recurrence relation). The results 
obtained by the proposed method are sufficiently accurate 
for practical applications, however, a modified definition 
to the proposed scheme is given m the later part of this 
chapter. The modified definition is formulated out of practice. 
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experience and this simplifies further the mathematical 
computations involved in the analysis of the proposed scheme. 
For convenience the following abbreviations are used for 
various schemes under study. 

Scheme 1 - proposed scheme ; 

Scheme 2 - Modified scheme (modified version of scheme 1) 

Scheme 3 - The method proposed in [101] , 

The results obtained by the Scheme 1,2 and 3 are compared 
with the exaot solution of the given system equation. The 
scheme 2 gives more accurate result for system whose solu- 
tion is known to be bounded in addition to a reduction in 
involved mathematical calculations. 

Example 2.7.1 : Linear (bounded/ unbounded) systc is : 

Consider the following linear difference equition with 
the perturbation parameter ji 

x(k+l) - x(k) + x(k-l) - \i x(lc-l) = 0 (2.45) 

(i) for 0 < \i < 1.0, tho solution decays to zero as k m *’ eo 

(ii) for ii < 0 and \i > 1.0, the response is unbounded, 
that is x(k) as k -**» 

(iii) for \i = 1.0 and for = 0.0, the solution is oscillatory 
(bounded) . 
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The application of perturbation methods for this type 
of snail perturbation problems has been found to give fairly 
good approximation to the exact solution. The eqn. (2.45) 
will now be studied for sma^l values of (i under unbounded 
as well as under damped conditions . 

The direct application of the discrete two time scaling 
approach given in eqns. (2.28) and (2.29) to (2.45) gives 

x(t)+1,t) - x(ti,t) + x(-n~l,T> -i- |i[2x(tj+^,t-4-} 

- 2x(n+£* ~ 2 x(ti- 4, *+£) + 2 x(r)-i^~t) 

- x(ri~l f t) ] + |i 2 { 2^ 1 x(t)+1,t) - (4^2+3) x(i),t) 

+ 2“^ x(t)-1,t) + 2x(ri,'r+l) + x(t},t;-- 1) + (2s 1 -2) 

[x(r)+-£-, t+4-) - z(t)H» *Hr) ~ x(r)-^, t+£) 

+ x(r}-4, t- 4-)]) + 0(p 3 ) = 0 . (2.46; 

Substituting eqn, (2.34) in eqn. (2.46) and collecting the 

coefficients of }i°, p 1 , p 2 , and eq.ua uing separately to 

zero, the following set of linear partial difference equations 
are obtained. 

* 

x 0 (t]+1, t) - x 0 (t),t) + x Q (r)-l, t) = 0 (2.47) 

x 1 (ti+1,'t) - x 1 ('n,v) + x 1 (ti-1,t) = -[2x 0 (r)+-^, T+-J-) 

- 2x 0 (r}+^., t- 4-) - 2x 0 ('-i- f t+^) + 2x 0 (r)-£, t-J-) - x 0 (t)-1,t)] 

( 2 . 48 ) 



63 


;; 2 (n+l, t) - x 2 (ri,'c) + x 2 (rj-l, t) = - [ 2x 1 (r}+£,x+-£-) 

- 2x 1 (ri+^ t T-4) - 2x 1 (‘n-4 f T+4-) + 2 Xi (ti-4-,t-4) 

- x 1 (t)-1,t)] - { 2 <*>iX 0 (Ti+l f T) - (4^+3) x 0 (t),t) 

+ 2* lXo (n-i,T) + 2x 0 (n,'t+l) + x 0 (r),T-l) + (2S-J-2) 

[x Q (riH» T+ i-) “ x 0 (tiH^-4-) “ x 0 (t)-4-,x+-Jt) 

+ 3c 0 (ir4»^-i)] > , (2.49) 

Bolutiun to eqn. (2.47) is 

x (tj,t) = A(x) cosOp + B(x) sinOr) (2.50) 

whore A(x) and B(x) are the slowly varying amplitude functions 
and 9 is the natural frequency of the system. Here 9 = w/3. 

Eqn. (2.50) is substituted on tlio right hand side of oqn. 
(2.48), using the relations given in oqn. (2.43), we get 

x^(rH*l,x) - x^f'Hjx) + x(r)-l,x) = -2[ i(x+^-) cos ^ 

+ B(x+£) sin cosOp +• 2[A(x+-|-) sin *r 

- B(t+|) cos §] Sj.n9ri + 2|_A(x~£-) cos ^ 

+ B(x~£) sin ■j] cos9r) - 2[A(x~4-) sin ^ 

- B(x— J-) cos sinOr) + 2|_A(x+-£) cos i| 

- B(x+£) sin cosQti + 2{_A(x-t-4 r ) sin v? 

+ B(x+-£-) cos sinOr) - 2[A(x-£) cos j| 

- B(x— £•) sm cosOri - 2[A(x~4r) sin § 

+ B(x~£-) cos sinOr) + [A(x) cos 9 - B(x) sin 9] 

cosOr) + (A(x) sin 9 + B(x) cos 9] sm9r). 
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to eliminate the secular terms on the right hand side of the 
above eqn. 

coefficient of cosGri = 0 
and coefficient of so lGr) = 0 

that is 

4B(t+|-) sm ^ - 4 B(t-£) sm £ = A(t) cos G - B(x) sin G 

4A(t~^) sm - 4 A(t~£-) sm ^ =-A(t) sin G - B(x) cos Q 

which can be rewritten as follows using the central difference 
operator 6, 

6B(t) = — — ~ — q |_A(t) cos 0 - B(x) sm 0] 

4 sm 

6A(t) = — “ — 2 [a(t) sm 0 + B(x) cos 0] (2,51) 

4 sm ^ 

The eqns. in (2.51) are first order linear coupled difference 
equations and these equations can bo solved easily as follows. 
With 0 = ic/3, the above equations can be rewritten as 

6B(x) = | [A(x) - f3 B(x) ] 

6A(v) =-i [V3 A(t ) + B(x) ] . (2.52) 

Differencing the first equation m (2.52) once 


6 2 B(t) = \ [ 6A ( x ) - V3 6B(x)j 
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Substituting for 6 B(x) from the second equation 
6 2 B(t) = 6 B(t) - ^ B(x) - H A(x). 

Again substituting for A(x) ’rom the first equation we obtain, 

6 2 B(t) + 6 B(t) + | B(t) = 0 . (2.53) 

As indicated in the remark at the end of the section 2.4, 
replacing 

6 2 B(t) = B(t+1) - 2B(t) + B(x-l) 
and 6 B(t) =$ 

the eqn. ( 2 . 53 ) can be rewritten as 


B(t+1) - ( 4 :^ 3 ) ®( T ) + ^"rf^ B(x-l) = 0 . 
The solution of (2.54) is of the form 
B(x) = R T [P cos cpr + Q sin 91 ;] 

where 



(2.54) 


(2.55) 


and P,Q are the constants to be determined from the initial 
conditions. 
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Then from the first expression m (2.52) 

A(t) = 4<5B(t) + V3 B(x). 

Substituting for B(t) and 6L(t) from the eq.ii. (2.55)» the 
following expression for A(t) is obtained. 

A(t) = R t [Q cos ft - P sin , (2.56) 

Combining the equations (2.56), (2.55) and (2.5C) we obtem 

x q (t),t) = R^[A cos (Or) - cpr) + B sin(0Ti - T'c)] 

Hence the approximate solution to the given equation upto the 
order \i is 

x(k) * x q (t},t) = [A cos Qjk + B sin O^k] (2.57) 

where 

0-^ = ( 0 - Jicp ) . 

Figures 2.2(a) and 2.2(b) give the comparison of the 
approximate solution given m (2.57) with the exact solution 
of the given equation for ji = 0.1 and |i = 0.3. It is observed 
that the approximate solution is closer to the exact solution 
for small \i, namely for \i = 0.1 for all tine and the deviation 
of the results increases as lc increases for larger \i m This 
situation is shown m Fig. 2.2(b). 






CXC- SO.J 

Proposed r 


r\j 

m 



b UNEAR JhDEE DAMPED SVSTEM (/J--0 3) 
x(k +1) -x (k 1 + x (k-l)-jux 'k-1 )-0 
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Figures 2. 3 (a.) and 2.3(b) show the unbounded response 
of a linear difference equation. Here again for snail |i 
the approximate result obtained by the proposed method is 
closer to the exact solution For large \i "n aopreciable 
diviation of the results is evident 

Now, the procedure to obtain the second correction term 
will be considered. Before proceeding with the calculations 
for second correction term it will be useful to say sonclhing 
about the initial conditions due to introduction of two inde- 
pendent time factors. 

let x(o) = x Q and x(l) = he given initial conditions 
fur the second order equation given m (2.45). Then from 
eqns. (2.27) and (2.34) we obtain the following initial con- 
ditions. 

(l) x Q (0,0) = x q 

x o (l,0) = 

(ei) x 1 (0,0) = 0 

x 1 (l,0) + x Q (0,l) - x Q (0,0) = 0 
(ni) x 2 (0,0) =0 

x 2 (l,0) + x 1 (0,l) - x 1 (0,0) + w 1 x {) (l,0) 

+ S 1 X 0 (0,1) - (% + Sj) x o (0,0) = 0 


and so on 




mss 


pip 


■ ,>! ' i 


mi ; : 


Him ' 


1 ill 


mM, 


MimW 







• * , ; ^v. 
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With x Q = 0.0 and x^ = 1.0, the generating solution for 
ji = 0.1, becomes 

x 0 (ti,t) =0 (0.63) X sm(1.0?7 k) (2.58) 

ulnre 

t = 0.1 k. 

Then the solution to eqn, (2.48) is given by 

x -L^ r ' ,T ) = C(t) cosOr) + D(t) sm©r] # (2,59) 

Substitution of eqns. (2.58) and (2.59) on the r j ght hold 
side of eqn. (2.49) leads to 

x 2 ('n + l»'c) - ^(t^t) + = “2lC(iH"£-) cos § D(T+£)oin 

cos©n + 2 [_C (t+ 4-) sm ^ - D(t+-£-) cos sinOr] 

+ 2[C(t— £-) cos ^ + D(t-~£-) sm cosOr) - 2 [_C(t— £-) sm £ 
- D(t-4-) cos sm0T) + 2[0 (t+-|-) coS § ~ D(t+-|-) sin §] 
cosOr) + 2[C(irt-Jr) sm + 0 ( t+ 4- ) cos sinOr) 

~ 2[0(t~J-) cos ^ - D(t— R sm *r] cosOr) - 2 [_C(t— ^) sin ^ 
+ D(t— : £-) cos sinOr) + [0 (t) cos 9 - I)(t) smO] cosOr) 

+ [a(t ) sm & + B(t) cos Qj smOr) - ? u coS 9 t 

T R % <r 

~ R /V"3 sinter] cosOr) - 2 w coS 9' t H smcpT] sinOr] 
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(4 w 2_ + 3 )y^ ^ siri V ) ' T coS 0'1 + (4 w j+ 3 ) 


2 T/ j-? v % 

yj R cosqDT smOr) + 2u-^iyy + R cos^t] 

r'C 

cos9r] - 2 W 3 _Ly5 go - R sirup tj smOrj + 


jj R x sin(<pT + 9) cos&ri - ^ R T cos^th-^) smQri + 

P 2 

rnrs R T sxn(ij)T -9) cos&n ~ • ■y"? p - R T cos (91; -9) sinQr] 


(2-2s 1 )l||^ R x sin (| - 9>t - |) - yfa R r sin(§-<pT+|) 


R t sin(|t- 9 T+|) — y^uR^ cin(^+ 9 T:-^) ]smQri + 

C0S (f “ 9^ - fr) - Y5 r rT cos (§ - 9 T + f) “ 

R 1 ' cos (| + 9 t + ^) + y|g- R t cos(^ + 9 T - |)] smUr). 


ELmirution of secular terns rosults in the following equations. 

o d 1 22, 0 cp 

2 w 1 4- y^ SlIKp - Sirup 4* (1 - ) [y^r sm(^ “ *r) 


“ tIr sin (f + f) + “H 


/ Q . cd \ 2 


2 oj 2 _ 4 u 2 . + ^ 2 R 1 

“Y3 — ( ^~Y3 ^ Y5 coe 9 7 jR G<_) ^ 9 (^* “ 

L~ Y3^ G0S (§ “ + y5r CuS (§ + 2 ^ + coS (§ + §) 

~ Y3R (2 ~ ^ = 0 . 

'.Those algebraic equations involving the two unknowns can now be 

Solved for w n s, and 

ly 1 
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a) 1 =-0.05 
s 1 = 1.247 • 

It is to be noted that the contribution, due to second corre- 
ction tern is more significant m ti c slow tine factor % than 
the fast tine factor r) since s 1 >> w , This con r radicts the 
statonent given by Morrison [104]- f <r coi bmuo j tine systen 

nanely the inclusion of the constants S do the definition 

n 

of slow tine scale being to account for t _c genorcl depondencc 
of the nonlinear function upon p. 

Table 2,1 shows the results of the linear discrete 
system described by oqn. (2.45) upto the order p . It is 
observed that the solution upto second correction terns is 
closer to the exact solution. This tabulated results are 
obtained for p = 0,1 in the eqn. (2.45). 

Table 2,1 

Solution of equation (2.45) with p = 0.1 

Discrete Exact Solution upto Solution upto 

time solution the order p the order p 2 

instant 
k 


1 

1.000 

1.000 

1.000 

2 

1.000 

0.995 

0.995 

3 

0.100 

0.078 

0.088 

4 

-0.800 

-0.829 

-0.808 


contd , » . * 
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Discrete 

bine 

instant k 

Exact 

solution 

Solution upto 
the order \i 

Solution upto 
the order \i 2 

5 

-0.890 

-0.896 

-0.883 

6 

-0.170 

0.135 

-0.150 

7 

0.631 

0.682 

0.64 5 

8 

0.784 

0.801 

0.7 7 8 

9 

0.216 

J.175 

0.192 

10 

-0.490 

0.556 

-0.511 

11 

-0.684 

- >.712 

-0,681 

12 

-0.243 

-0. 202 

- 0 , 217 

13 

0.372 

0. 449 

0,398 

14 

0.591 

0.63C 

0.591 

15 

0.256 

0. 218 

0,229 

16 

-0.276 

-0. >58 

-0.804 

17 

-0.506 

-0.655 

-0.509 

18 

-0.258 

-0. 225 

- 0. 732 

19 

0.198 

0.281 

0.228 

20 

0.430 

0.485 

0.436 

21 

0.252 

0.226 

0. 228 

22 

-0.135 

- 0.217 

-0.166 

23 

-0.362 

-0.422 

-0.370 

24 

-0.240 

0.222 

-0.219 

25 

0.085 

0.164 

0.116 

26 

0.302 

0.365 

0.312 

27 

0.225 

0.214 

0.207 

28 

-0.046 

0.120 

-0.076 

29 

-0.249 

-0.315 

-0.262 

30 

-0.207 

-0.203 

-0,192 

31 

0.017 

0.085 

0.045 

32 

0.203 

0.269 

0.218 

33 

0.188 

0.191 

0.176 

34 

0.005 

-0.056 

-(,.021 

35 

-0.169 

-0.229 

-0.179 
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In the following section a linear oscillatory type of 
system with small perturbation is considered and the analysis 
is carried out with the proposed scheme and the results are 
verified by computer simulation. 

Example 2.7.2 : linear discrete oscillator : 

Cunsidcr the linear discrete oscillator described by : 

x(k+l) - x(k) + x(k-l) + (3 x(k) = 0 (2.60) 

Depending on the numerical value of (3, this equation can 
be rewritten as shown below . 

(l) for valves -1.0 < p < 1.0, (3 can bo taken as a 

small parameter p, that is p = p and the equation (2.60) 
assumes the form 

x(k+l) - x(k) + x(k-l) + p x(k) = 0 (2.61) 

(ii) for values 1 < p < 3, p can be rewritten as 
(j = a + p, so that p = p-oc is always less than unity. Under 
this condition the system equation can be rewritten as 

x(k-i-l) - (1-a) x(k) + x(k-l) + p x(k) = 0 . (2.62) 

how the eqns. (2.61) and (2.62) are of the form to apply 
two variable expansion scheme. Considering tho eqn. (2.61), 
the basic Solution is 

x 0 (r),x) =A(x) cos Op + B(t) sinQr) (2.63) 


with 0 = te/3. 
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file right hand side of difference equation relating 
(first correction tern) is 

= 2 x 0 (r}-»4, t+^) - 2x o (ri+4-, T “i) ~ 2 ^ 0 (r]-4, x+±) 

+ 2x 0 (n-i, ) + x 0 (t},t). (2.64) 

Substituting (2.63) m (2.64) and equating the coefficients 
uf cusQq and sinGq separately to zero, the following equations 
are obtained . 

SB(t) = - 1 — g A (t ) (2.65) 

4 sm ^ 

6A(t) = — — g- B(t) (2.66) 

4 sm ^ 

Differencing eqn. (2.65) once on both sides with 6 

<5 2 B(t) = 6k(x) 

4 sm | 

Ooiibining the eqn. (2.66) with the above equation, we got 

6 2 B ( t ) + 0. 

16 sm ^ 

With 9 = %/3, and expanding, the following linear difference 
equation is obtained 

B(t+1) - 1.75 B(t) + B(r-l) = 0. (2.67) 


The solution to the above equation is of the forn 
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B(t) = p cosfpx + q sincpT 


( 2 . 68 ) 


whore 




(1.75/2) 


) 


and p and q are the constants to be evaluated knowing the 
initial conditions. Then iron eqn. (2.65) 


a(t) = p sirup t - q cus <pr . (2.69) 

from equations (2.68), (2.69) and (2.63), the approximate 
solution upto the order p is given by 

^(rjj'O = A cos Qjk + B sm O^k (2.70) 

whore 9^ = 0 + <pp. and A and B are the constants. 

The solution given m oqn. (2.70) is compared with the 
exact solution of the given equation m figures 2.4(a) to 
2.4(d) for various values of p. 

It is observed that the phase difference betwoen the 
solutions increases as p increases, but the magnitude of the 
solution remains unaltered. 

Example 2.7.3 : Weakly nonlinear difference equation 

In this example a discrete time system described by a 
weakly nonlinear difference equation is considered. The 
discrete version of the Buffing equation with weak forcing 





mm 




ifp 




v* >%?!' 


mm 


';: "4-w 


i§g4f v 

|l|§i 11 ' i 


-V ^ 


• - ;: : 


■ .■■••:•-■ :■ i‘ 

' *. 

S-sfcVi.: 


; -r! ;• ';;K: 


; %J ‘; - : 4 f 











: V" ; 


I g§|SJ 
Mil I 


« i 


v . 


.;; f i-.-i i v- ■■ ^ 

Mm H 


’ 1 ‘ 1 

'■'■’JS'figi 






• •'" I II 


<s $ Jit | & 











83 


Is under investigation. Vote that. this enables one to nake 
a comparison of the discrete model analysis with the continuous 
time system solutions. The example (a) discusses the situation 
where the continuous time version is dissipation free and the 
example (b) treats the oase where there is weak vlsoous damp- 
ing. The discrete models are derived in Appendix 3 and the 
analysis of these models are carried out here. 

(a) Duffing Equation (Dissipation free): 

The discrete model is 

x(k+l) + Xx(k) + x(k-l) + p, y x^(k) = V& oos u> k , (2.71) 

Adding and subtracting 2oos u x(k) in the above equation and 
rewriting . 

x(Ih-X) - 2 008 » x(k) + x(k-l) - coso k - y x 3 (k) - Kt(k)] 

(2.72) 

where M is the detuning parameter and is defined as 
l&M ^ X + 2 oos u • 

It is to be noted that the addition of the term 2oos « x(k) 
in the above equation assures bounded linear response. 

The application of the two time perturbation prooedure to 
eqn. (2.72) leads to the following set of linear difference 
equations 

x 0 (tj+l*v) - 2 oos « x q (t|,t) + x 0 (rj-l t T) ■ 0 


(2.73) 
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x 1 (tj+1,t) - 2 oos u x 1 (t),'c) + x^tj-J-.t) = P oos« rj - y*^(i»**0 

- Mx 0 (t),t) - 2[x 0 (Ti+4r»*rt4 - 

- x 0 (t)-4 f T-4.) + . (2.74) 

The solution x 0 (t),t) obtained fron equation (2.73) is . 

x 0 (t),t) * A(t) oos" tj + B(t) sinw tj 

Substitution of the generating solution using (2.43) in the 
right hand bant side of oqn. (2.74) yields 

X^(r)+l f x) - 20080 ^(TJ,T) + X^TJ-ljX) = Poosw tj 

- y[A(t) oos «t) + B(t) sinw ti]^ - M[A(?) cosw tj 

+ B(t) sin« t}] - 2[A(t+^) oos ^ + B(t+ 40 sin^]ooew tj 

- 2[A(‘d4>) sin ? tj - B(x+i-) oos ig] sin" tj + 2[A(t-£) . 

oosw/2 + B(-c-±) sin" /2] cos® tj - 2[A(t- 4) Bin«/2 

- B(x~£) cos®/2] sin" rj + 2| > A(t+J.) oos®/2 

- B(t+4-) sinw /2] ooswt) + 2[A(t+^) sin®/2 

+ B(t+£’) cos ®/2] sinUT) - 2[A(t-£) oos w/2 

- B(x-fc) sinu /2] oosw tj - 2[A(t— $■) sin w/2 

+ B(t-J-) oos w/ 2] sinwt) (2.75) 
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In the above, nonlinear term oan be rewritten as 

[A(x) oos® t) + B(x) sins rj]^ - A^(x) + J A(t)B 2 (t)]oos«9 tj 

+ [f ^(t) + l A 2 (t) B(t)] slnu i) 

+ Non-secular terns 9 (2*76) 

Substituting the eqn. (2.76) In eqn. (2,75) and collecting the 
ooofficients of oosut) and sin u ri and equating then separately to 
zero the following equations for A(t) and B(x) are obtained ; 

«a(t) - [| 7 »(■») U 2 (t) + b 2 (t))+ m b(-c)J 

“ (t) “ X a ' lniTTS £* - | y a (t) U 2 (*> + b 2 W)- M A(t)J 

(2.77) 

where 6 A(t) = A(x + £) - A(t- 40 and likewise 6 B(t). 

If a steady state solution for a(t) and B(x) exists then 

6A(?) = 6B(x) a o 
and A(t) = A, B(x) = B 

where A and B ore the steady state values of A(x) and B(t) 
respectively. 

Then under steady state condition 

B[| 7 (A 2 + B 2 ) + M] « 0 (2.78) 

F - A y (A 2 + B 2 ) + M] - 0 • 


(2.79) 
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From eqn. (2.78) 

either B *= 0 

or M = - | y (A 2 + B 2 ) . 

Ouiibining the second condition with eqn* (2.79) » we get 

F * 0, corresponds to force free situation* 

So ruling out this condition, the first condition namely 
B = 0, with eqn. (2.79) gives 

F * | y A 5 + Mi ■ 

from which M«J-^yA 2 . 

Substitution of H yields 

cos “ * |i/2 (f - | Y A 2 ) - X /2 . (2.80) 

Eqn. (2.80) shows the output amplitude variation with the 
irqput frequency. 

Case 1 s If y ■ 0*0 (linear spring) 
cos “ ■ - X/2 

gives the linear system response and is plotted in Fig. 2.5 
for various values of the input amplitude F. 

Case 2 t If y = 1. 0 (Hard spring) 


cos u 
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9he response characteristics under this situation is shown 
in Fig. 2.6 for various F. There is a striking resemblance 
to the familiar response curves of the continuous Duffing 
equation without damping [5]. 

Qubo 5 t If y ■ -1*0 (Soft spring) 

ooe m - | [| + \ A 2 ] - * /2. 

The variation of the amplitude with input frequency for 
various values of input amplitude is plotted in Fig. 2.7. 

In all the three oases the response curve corresponding 
to F * o shown by dotted lines is known as * backbone* . It 
is interesting to observe that the familiar bending 

of the response curves for y £ 0 about the backbone 

arises for this discrete time system also* 

(b) Duffing equations (with weak damping) 

For the case where there is weak viscous damping the 
discrete model arrived at is 

x(k ) + X x(k) + x(k-l) + (i C[x(k+£) - x(k-£)] + * 

11 ^ cos a k - F 2 sin « k] 9 (2.81) 

where C is the damping coefficient assumed positive. 

Adding and subtracting 2 oos u ac(k) on both sides of 
eqn. (2.81) we obtain 






■■■■■■■I 


iMIbm. 


iSSiill lilllll 

;M : :.;:jaj|f 


MBa ilfM llj 


*■■■:■ : 


/ v. ■ , \ ■ ■/ ■■ 

*1 • - ; m::i 

. **' ' ‘ > Ws**^ t*, * ’ »%, ?*■' r* .*» <s 

a- • ■»;* 44 ?' ■■■., '.«-■ V ' ‘;,» . i . •:*. 

■ r *+« * ., '; - ,, * N ,. ■• 

' ’" . -* >****'•** •/ . fr\ ' . 4 

j ** . *\ _AS *, ,*"■ *_ T ' .. , •, 


;$ht- 


'«4 j,n o r ^? 

.• . T? 

v . . . - tL 




; , if 

: .: .- • . . 

■■ -'ya. 

* 

; : ■ 








v'VV? J ' 

'■ 


1 - 


>•• - •' 






91 


x(k+l) - 2 oob • x(k) + x(k-l) + |i M x(k) + |i 0 [x(k*£)- x(k— J-] 

+ |i Y oob « k - Pg sin“ k • (2.82) 

Proceeding as before the basic solution x q (t| v t) is again ob- 
tained as 

x q (t) 9 t) = A(t) costa t) + B(x ) sins tj • 

Examining the terns of order |i leads to the following 
variation for A(t) and B(t) aftor supression of secular terms 

6A(x) “ Vain « n t*2 + H B(t) ’ 2 0 A ( x ) aln “ / 2 + } YB(x) 

(A 2 (x) + B 2 (x)] 

8 B(t) . 4 -^; a y ; - K A(x) - 2 0 B(x) ai*®/2 - | r A(x) 

(A 2 (x) + B 2 (x»] 

from which the steady state behaviour determines 
P x = Hi + 2 0 B sin « /2 + | y A (A 2 + B 2 ) 

P 2 « -MB + 2 0 A slnu/2 - | y B(A 2 + B 2 ). (2.83) 

lot 

F 2 * F 2 + p| 

where F is the magnitude of the forcing function, that is 

F 2 » [MA + 2 0 B Bin « /2 + | y 1 (A 2 + B 2 )] 2 

+ [MB - 2 0 A sin M /2 + JyB(A 2 + B 2 )] 2 
= R 2 [M 2 + 40 2 Sin 2 «/2 + fa yhfi + \ MyR 2 ] 


(2.84) 
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With 0*0 and F^ * 0, the eqn. (2.83) reduces to 

r x - » + | ya(a 2 + b 2 ) 

0 - MB + | fB(A 2 + B Z) 

which are same as given in eqns. ( 2.78) and (2.79) . 

Due to the presence of the weak damping term, closed 
response curves for the Duffing equation are shown in Figs*. 
2*8 to 2*10. It is observed from the figures, for a given 
amplitude F, each response curve has a vertical tangent as 
well as horizontal tangent. She equations of the curves des- 
cribing the locus of the vertical and horizontal ta agencies 
are derived as follows* 

Locus of the horizontal tangency points : 

Consider the response chara ct eristics given in (2.84), 
that is 

F 2 * R 2 [M 2 + 4^ sin 2 " J2 + £5 + \ MyR 2 ] . 

Here the variables w and R are assumed to be differentiable 
with respect to one another. 

She locus of the vortical tangency is obtained by 
Imposing the well known condition 

& -0. 

For convenience the eqn. (2.84) can be rewritten as 
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^2 * (M + ^ yE 2 ) 2 + 4<?sin 2 »/2 . 

Differentiating with respect to u , (noting that the amplitude 
E is a function u ) and setting g— « 0, ve obtain 


0 

where M 


2[m + | y a 2 ] Sjf + 2< ^ Bia “ 

2 oobihX 


and « -2/n sin a . 


Substituting in the above equation, we get 


008(0 



( 2 . 86 ) 


Fig. 2.11 shows a plot of the horizontal loci for the cases 
C = 0.0 and C = 1.0 9 that is without and with weak damping. 

Locus of the vertioal tangency points : 

The locus of the vertioal tangents oan also be easily 
obtained from the eqn, (2.84). 

The condition to be satisfied here is 

Sft-o- 

Differentiating the eqn. (2.84) with respect to E. and Inserting 
the above condition, we get 

(M + \ yH 2 ) (M + J y H ? ) + 4 C?sin 2 « /2 - 0 . 
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Substituting for H and simplifying, the following equation, 
representing the locus of the vertical tangents is obtained • 

oos 2 « - ^ oosot + ^ * 0 (2.87) 

■ where 

■ ■¥ - 1 ^ - x 

^ J [2|» 2 o 2 + + 3 |iy * R 2 +' * 2 J . 


Fig. 2.12 shows a plot of the vertical tangents looi for the 
cases C = 0.0 and' C = 1.0. It is to be noted that for small 
values of 0 the vertical tangents looi are close to the curves 
defined by 

H + | y R 2 - 0 
M + ^ y R 2 a 0. 


Jump Phenomenon : 

Figure 2.13 indicates the amplitude -frequency response 
characteristics with weak damping for a particular input 
amplitude F when y = 1.0. The 'jump phenomenon' that occurs 
for variable Inputs frequency at a constant amplitude is also 
shown in the same figure. For increasing « , 'jump' occurs 
from a to b and for decreasing a , 'jump' occurs from o to 
d. These properties are well known for continuous time systems 
and the above analysis provides the discrete time counterpart 
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as veil as supporting the validity of the discrete multiple 
time perturbation analysis technique. The stability analysis 
of the output amplitude function can be carried out by locating 
the singular points In eqn. (2.83) Tor specified Input ampli- 
tude and frequency. The above stability study Is straight 
forward and the results are expected to be of the same form 
as In continuous time systems and hence not considered In 
the present analysis. 

Example 2.7.4 Von dor Pol Type Weakly Nonlinear Equation : 

A Tan dor Pol type nonlinear discrete oe dilator is con- 
sidered and Its limit cycle behaviour is studied. The system 
is described by the following difference equation, 

x(k+l) + x(k-l) + |i[l - x 2 (k-l)] x(k) - 0. (2.88) 

The base solution is 

x 0 (t),t) ■ A(t) cos On + B(t) sin On 
with 0 m 1 t/ 2 . 

Considering the solution of x^(n*'0 » the following equa- 
tions for amplitude functions A(t) and B(v) are obtained . 

2f2 SB(x) - A(x) [1 - B 2 (x)] 

2V2 61(x) - -B(x) [1 - A 2 (x)J . (2.89) 

Under steady state condition 
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A(l - B 2 ) * 0 
B(1 - A 2 ) - 0 

from which A « 1, B = 1 are the steady state amplitude* 

This Implies that the Vondor Pol type system described 
in eqn* (2.88) exhibits a limit cycle oscillation of amplitude 
unity. The stability of such limit cycles are analysed through 
variational technique as follows : 

Before studying the limit cycle stability, it is useful 
to study the singular points of the system. The real sin gular 
point is easily determined to be the origin x = 0* The sta- 
bility of this singular point is now analysed by considering a 
small variation about it* let the variation be x « x* + u, 
where x* is the singular point* Substituting in the given 
system the following variational equation is obtained, 

u(kt-l) + u(k-l) - jiu(k) « 0 . 

The solution to the above equation is oscillatory, which 
implies that the singular point is a center type of singularity* 
This shows the given system exhibits oscillations for the 
initial conditions which are very close to the origin* 

The stability of the limit cycle oscillation of unit 
magnitude is studied by considering a small perturbation 
around unity, namely 
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a(t) = i + a a(t) 
b(t) « i + a b(*0 . 

Substituting these sanations in eon. (2.89) and neglecting 
higher order terms the following variational equations are 
obtained i 

+ 3 B(x+1) ? 3B(x) + - 3b(t) * 0 

“ Y2 

+ &A(t+1) +3 A(t) + i-8A(t) «0, 

V2 

It is evident from the above equations aB(t) ♦ o as t 
and 3 A(x)-*® as % **• , which shows that the limit cycle is 
of semistable type. This has been verified by simulating the 
given system equation. The figure given belows illustrates 
the stability aspects of limit cycle oscillation in the incre- 
mental phase plane. In the figure, A is the initial condition 


Z(Mi) 



A r 

“1- 



L 

J 

XCht) 





very close to the origin, B and C are the initial conditions 
inside and outside the periodio limit cycle os dilation (here 
period 4) • 
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2.8 Modified Definition of Two lime Scaling Method s 

As mentioned earlier the proposed scheme works veil 
for all types of equations In a Class of discrete time systems. 
This has been shown In the previous section by considering 
various kinds of linear as veil as nonlinear equations. How- 
ever, a modified definition over the initially proposed scheme 
Is if* Introduced here which, has a smaller number of mathematical 
calculations and obtains solutions very dose to the exact 
solutions for certain equations whose solutions are known to be 
bounded. 

The two independent time scales, namely T), the fast time 
scale and x, the slov time scale are defined as given In eqn. 
(2.19) whereas the dependent variable x(k) Is expanded In tosas 
of t),x as shown belov : 

x(k+l) * x(t)+1,t) + h[x(t)+1 # t+£) - x(rj+l, x-$)] 

x(k-l) a x(t)-1,x) - |i[x(T]-l f T+£) - x(tj-1, x-£)] (2.90) 

The expressions in eqn. (2.90) give the modified definition 
for the proposed scheme. This modified scheme (scheme 2) is 
applied, to analyse linear as well as a class of nonlinear dis- 
crete time systems and the details of the subsequent analysis 
are exactly similar to those for the proposed scheme (scheme 1). 
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Hie results obtained by the schemes 1,2 and 3 are verified . 
with the exact solution of the system equation and a compara- 
tive study of these schemes is given in the next section. 

2.9 Comparison : 

The scheme 2, given in the above section and the sohsme 
3 » proposed in [101] are applied to the examples 2.7.1 and 
2.7. 2 and the results are compared with the saherne 1 and the 
exact solution obtained through simulation of the system equa- 
tion. The plots showing the above comparison are shown in 
Figs. 2.14 to 2.16. As mentioned earlier the scheme 2 is better 
than the other schemes when the solution to the system is 
known to be bounded. Hie steady state equations obtained for 
nonlinear systems, by application of all these schemes, are the 
same and this has been verified by adopting the above schemes 
to the examples 2.7.3 and 2.7.4. 

2.10 Conclusion : 

The multiple scale perturbatlunal approach which has 
been used so far for a class of nonlinear differential and 
partial differential equations has been successfully employed 
for analysis of a class of nonlinear discrete time systems 
described by difference equations. The discrete multiple 
scale perturbatlonal scheme has been obtained in a systematic 
way using the known properties of the finite difference ope- 
rators. and nonlinear examples supporting the usefulness 
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of the proposed method have been presented. A modification 
over the proposed scheme is suggested. This modified scheme 
involves lesser mathematical computations and is more suitable 
when the solution to the given nonlinear difference equation 
is bounded. A comparison of these methods with the scheme 3 
as veil as with the exact solution of the system equation has 
also been presented. 



CHAPTER. 3 


STBOHOU FORCE]) NONLH9SAR MEFERBTCB EQOATIQHS 
3.1 Introduction t 

Weakly nonlinear discrete time system under free and 
weakly forced situations were considered In the previous 
chapter and several significant results were obtained. In 
this chapter the study is concentrated on weakly nonlinear 
discrete time system with strong external input. The study 
of second order continuous time system under strong forcing 
has been carried out by many persons. When the forcing is 

i 

strong, various additional phenomena occur that are not 
present in weakly farced or in force free systems, specifi- 
cally, the generation of super/ subharmonio oscillations and 
frequency entrainment. It is well known that frequency en- 
trainment is possible in a self exited system with external 
excitation [4*5]. 

Hayashi [4] has investigated subharmonlc solutions for 
polynomial nonlinearity by the method of harmonic balance. 
The resulting algebraic equations were solved by Iterative 
method to obtain the required solution. Kronauer and Musa 
[23] have also reported some general expressions for all 
subharmonios using multiple time perturbation approach when 
single input is applied. Tiwari and Sibramanlan [22] have 
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reoently obtained some necessary conditions, for 'various 
super, sub and ultxasubhaxmonic synchronizations in a class 
of seoond order differential equation with multiple excita- 
tions. iomas [114] Investigated the existence of 3/2 ultra- 
subhaxmonlos In a Puffing system through the method of haxmonio 
balance. Stanlsio et al [113] presented a general perturba- 
tion method for the study of sub and super harmonic solutions 
of weakly nonlinear second order differential equation without 
damping. Whereas recently Riganti [116] has studied the sub- 
harmonic solutions of weakly damped Duffing equation with 
large nonlinearity. 

In this chapter, the method proposed in the previous 
chapter for analysing the force free and weakly forced seoond 
order difference equation, has been adapted to investigate 
additional phenomena that are generally present in strongly 
forced situation. Here again the nonlinearity considered for 
the study is of polynomial type. A Puffing type nonlinear 
discrete osoillator is considered and various possible super 
and subharmonio oscillations due to a strong external Input 
are worked out. The well known technique of Harmonic Balance 
has also been employed to study the above phenomena. £he 

t 

I results deduoed from the proposed methods are compared with 
the ones obtained by computer simulation of the system 
equations. 
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3*2 System Description 

Consider a system which is described by the following 
difference equation 

x(ta-l) + ocx(k) + px(k-l) -i- pf (x(k), x(k-l) ) » ? cos » k, 

(3.1) 

where 

a and p are the system parameters, p is the perturbation, 
parameter (generally small) and F, u are the amplitude and 
frequency of the input forcing function. The system in eqn. 

(3.1) is said to be strongly forced if the amplitude P is 
of 0(1). 

It is to be noted that, for the location of [a,p] in the 
stability triangle, the linear system response is asymptotically 
stable/ where as the linear response is periodic for [a,P] 
values on the side AC (except at points A and C) of the stability 
triangle ACC shown in Appendix A. Vithout loss of generality 
the value of p can be taken as unity for bounded periodic res- 
ponse with | a | < 2.0. A weak damping, to the system des- 
cribed in eqn. (3.1) can be introduced by assuming p nearly 
equal to unity. With this description the analysis of exi- 
stence of possible super/sub harmonic solutions is given below. 
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3.3 Analysis : 

Adding and subtracting x(k-l) in eqn. (3.1) , It oan be 
rewritten as 

x(k+l) + ocx(k) + x(k-l) + pNx(k^l) + pf(x(k), x(k-l) 

= F cos m k f (3.2) 

where 

liflf s (p-1), and R is a damping detuning parameter. It 
is to be observed that for p = 1, the detuning parameter 
N « 0 and the given system is free from damping. 

Due to forcing, the linear solution, for p as o is obtained 
as follows s 

She linear foroed system takes the form 

x(k+l) + ox(k) + x(k-l) as F cosu k. (3.3) 

Let the steady state solution be 

x(k) «= S oosuk + Q sino t, (3*4) 

where £ and Q are the constants to be determined from the input 
periodic function. 

Then, 

x(k+l) as (p oos e + Q sins ) cosu k - '■* 

(P sins - Q cosu ) sinu k 
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x(k-l) = (P oos « - Q sin" ) oos ® k + 

(P sin ® + Q oos « ) sin® k, 

Substituting these relations in eqn. (5*5) and simplifying 
P(o + 2 oos “ ) ■ P 

Q(a + 2 cos ®) » 0, (3.5) 

the second expression in eqn. (3.5) implies 
Q as o and (a + 2 oos u ) 0. 

Then the linear forced response is given by 

x(k) a P cos® k. 

This response for the linear system suggests the following 
transformation for the complete nonlinear system (3.2) 

x(k) ■ y(k) + P oos ® k, (3.6) 

from which 

x(k+l) s y(k+l) + P oos® cos® k - P sin® • sin® k 

x(k-l) a y(k-l) + P cos® .oos® k + P sin® • sin® k # 

Substitution of these expressions in (3.2) results 

y(k+l) + cty(k) + y(k-l) + iiHy(to-l) + |if(y(k) + 

+ P oos M k, y(k-l) + P oos w . oos w k 

+ P sin u sin ® k) a-jiN (P cos ® • oos® k 

+ P sin “ • sin*® k). 


(3.7) 
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The eqn. (3.7) Is a weakly nonlinear discrete time system 
with weak forcing. Adding and subtracting 2 cos ~ u ,y(k) in 
eqn. (3*7) and simplifying, 

y(k+l) - 2 cos ,y(k) + y(k-l) + wy(k) + pHSr(k-l) 

+ |tf(y(k) + P cos « k, y(k-l) 4- P cos u . oosu k 

4* P sin u . sin uk) > ~|iB (P oos u . oosu k 

4- P sin u sin to k), (3*8) 

where Y is the frequency detuning parameter defined as 

YU ■ 2 cos - • + a, (3.9) 

and n» m are positive integers. 

The solution to eqn. (3.8) oan he obtained by applying 
the discrete multiple time perturbation technique proposed 
in the previous chapter for a class of discrete time systems. 

Ve shall use this technique and will represent the solution of 
(3.8) by an asymptotic expansion involving two time variables : 

y(n» u) =* y 0 (r» t t) 4- p y^fa,*) 4- liS^Olf*) + ••• (3.10) 

where t) is the fast time and v is the slow time and are treated 
as independent. The definition of these independent scales 
and the expansion of the variables x(k), x(k-fl) and x(k-l) in 
terms of these independent variables are given in the previous 
chapter. It is to be noted that in the first order approxima- 
tion to the solution the Independent scales assume the 



following simple form 


U7 


n -k 

X m Ilk. 

Then 

yOO * y(ti»T) 

y(k+l) + y(k-l) = y(n+l,x) + y(Tj-l,x) + 2|i[y(iH£» < c+£) ~ 

y(n+4-»‘t-i) - yto-fr,***) + y(ri-M-4)] . (3.n) 

Substituting the eqn. (3*11) in eqn. (3.8) 
y(T|+l,T) - 2 oob gu y(ti,x) + y(n-l,T) + 2ii[y(nr*4,x+£) 

- y(tiH» T *4) - yCTH-*^) + yCtH^-i)] + \tyj(n 9 %) 

+ iiNyfa-l,?) + |if[y(rj,T) + P oos» tj, y(t|-l f T) 

+ u( •) + p 008 M 008 tt T) + P slnuslno rj] a 

- pH(P 008 ». ooso) i) + P sin u • sinoi tj). (3*12) 

Combining the eqns. (3.10) and (3.12) and oollecting 
terms of like order in p, the following equations are obtained 

|i° terms x y 0 (rH-l # x) - 2 ooe § a y 0 (r).*0 + y 0 (tj-l,T) - 0 

(3.13) 

|i terms x y^rn-l,*) - 2 oos ga y 1 (Ti,x) + y^ti-l,*) - 
-N(P coso> . oos«n + P sin siitfo tj) 

- f(y 0 (v)»*) + P oosut), y 0 (T>-l,x) + P eosa . oosa n 
+ P sina . sina rj) - N y 0 (r)-l t T) - Try 0 (T),x) 

- 2[y 0 (r)+*#xH) - y 0 Ol+Jr» T “i) - y o 0Hfc»'*H) 

+ y 0 OHr» xwj)] # 


(3.U) 
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She solution to eqn. (3*13) 1 b known as basio or generating 
solution and Is given lay 

y 0 (T|**0 - A(t) oos §»t| + B(*0 sin jjaTN (3.15) 

Hote that the response 7 0 (t) v t) Is purely oscillatory and this 
arises primarily because of the particular from of w that 
Is defined In ( 5 « 9 ). 

Where 

A and B represent the amplitude functions of 

I) SubkarmonioB, if n/xa < 1.0 for n * 1 and m > 2 

II) Superharmonio , If n/m > 1.0 for xn = 1 and n > 2 and 

ill) Ultra subharmoni ob , if n/m < 1.0 for n and xa relative 

> 

prime Integers. 

The solution in (3.15) is then substituted into the 
right hand side of eqn. ( 3 . 14 ) and the secular terms are 
suppressed as follows. 

Let = § “ , then 

y 0 ( T l+i»‘ M 4 ) • [*(*+£) cos + B(-c-Hf) sin -J] cos u-jTi 

W*l ® 1 

- [A(vt£) sin -g - B(v+£) oos HgJ sin®^ 
y 0 ( T )~i»‘ p Hr) * C A (‘*H) oos - B(xH) sin - 5 ] eosw^ 

+ [A(t+£) sin^Jj + B(v*£) oos ^3 sin 
the right hand side of eqn. (3.14) becomes 
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BBS s -N(P oos cs ooso rj + P 8 into , slaw n) - oos a^r) 

- sin -E[A(t) cos®j~ B(t) sinu oos » ~t] 

- H[A(t) sin \ + B(t) oos “J sin®^ - y A(t) oos®^ 

Y B(x) sin ttjt) - 2 [A(t+£) cos ttl /2 + B(t+^) sln“ /2] 
oos + 2 [A(t+£) Sin - B(*pH^) oos sin « jr\ 

+ 2[A(t-£) oos -g + 3(x~i) sin OOS 8 jTJ 

- 2[A(t- 4) 8 in "4 - B(r-i) 008*^3 sin® jtj 

+ 2[j(t+ 4) oos - B(?t£) sin*^] oos®jT| + 2[A(-p*$) 
sin + B(t+£) oos -^] ela®^ - 2 [A(t-£) oos -ij 

Aft 

- B(t-&) sin -^J cos jT| - 2 [A(t- 4) sin + B(t-£) cos 
sin ®jTj« 

Equating tbe coefficients of oos ®^tj and sin ® ^tj te rms separa— 
tely to zero (secular term elimination) 

oob « jT) terms s 

-a^ - I(A(t) oos oi ^ - B(t) sin® ^) - yA(t) - 46B(t) 

• ■ CO-j 

sin ‘*|*0 


that is. 


“i 

4 sin -g 6 B(t) 


-a^ - 1 T(A(t) oosm^ - B(t) sin u^-yA^) 

( 5 . 16 ) 


sin w jt) terms i 
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-b x - N(A(t) sin*^ + B(x) ooa« 1 ) - yB(t) + 46A(t) sin^ = 0 
that Is 

®1 

^ sin 6 A(t) « b^ + IT (A(t) slnu^ + B(t) oos « , ) 

+ yB(t). (3.17) 

Bqas. (3.16) and (3*17) are first order coupled difference 
equations (usually nonlinear) In A(t) and B(t). In eqns. 
(3.16) and (3*17) » a^ and b^ are the coefficients of the 
fund a mental components cos u ^ri and Bln » ^r) of the w«ni in ear 
function f in eqn. (3*14) and 6 is the central difference 
operator defined as 

61(t) = A(t + i) - A(t - ±) 

and likewise 6 B(t). 

Under steady state condition 

6A(t) = 6(B) (t) ■ 0 and A(x) = A, B(t) = B 

where A and B are the steady state -values of A(t) and B(t) 
respectively. Then the equations (3*16) and (3*17) take Idle 
form 

a^ + H(A cos " B sin “|) + yA ■ 0 

b^ + Iff (A sine j + B oos w ,) + yB * 0 (3.18) 

Eqn. (3.18) gives the general steady state condition on the 
amplitude functions A(t) and B(t) for a general nonlinear 
function f • 
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la the next section a particular form for the nonlinear 
function f is assumed and the possible super/sub harmonic 
solutions are investigated* 

3.4 Example i 

Consider the following nonlinear difference equation 
with weak oubio jfaonlinearity . This discrete oscillator 
may be considered as a Duffing type with weak damping. 

x(k+l) - x(k) + px(k-l) + nx^(k) * F cos u k (3.19) 

Comparing (3.19) with (3.1) 

a * -1.0 

*( *)- k). 

It is to be noted that as stated earlier in this chapter, 

P may be taken to be slightly less than unity to introduce 
small damping in the system. 

The nonlinearity given in (3.19) takes the following 
form due to the transformation of variable [eqn. (3.6)]. 

x’Ot) - [jr(k) + P 008“ k]^ 

■ y^(k) + 3Py 2 (k) oosu k + 3P ^yCk) oos 2 s k 

+ pW ® k. 

Due to the introduotion of two Independent time scales the 
above expression can be rewritten in the following form 
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- y 3 (t|,T) + 3ly 2 (t|,T) cos w t) + 3P 2 y(rj f T)oos 2 » ti 
+ P 3 00S 3 e T). 

Then using eon. (3. 10) the above eapression oan be rewritten 
as 

a^O)***) * y^rwO + 3Py 2 (T),T) cos u t) + 3P^T 0 (ii,t) oos 2 e tj 

+ P 3 oos 3 e tj + 0(ji). 

Substituting for y 0 (T],x) from eqn. (3.15) 

y^(ti,T) + SP^Ol,*) cose t) + 3sV 0 (n.t) cos 2 ## TJ + P 3 OOS 3 e TJ 

* [A(t) oo se jTj + B(t) sine^] 3 + 3P cose cos « n 

+ B(t) sin + 3P 2 cos 2 erj[A(T) cos e^ 

+ B(t) sin e^Tj] + P 3 cosset). (3.20) 

BHS of (3.20) * A 3 (t)cob^ e^n + 3A 2 (t) B(t) oos 2 » sinB^T) 

+ 3A(t) B 2 (t) oos sin 2 + B 3 (t) sin 3 ® ^ 

2 2 

+ 3P oose t)[A (t) oos + 2A.(t) B(t) oose^r). 
sin b^ti + B 2 (t) sin 2 b^tj] + 3P 2 oos^» t)[A(t) oos 
+ B(t) sin e^tj] + P 3 eos 3 o n. 

Using well known trigonometric identities, the above expre- 
ssion can be rewritten as 
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■ + |a(x)B 2 (x) + |p 2 A(t)J ooBttjtj + [|b 5 (t) 

+ Ja 2 (t)B(t) + %?hix)} Bln ttjTJ + [J«L 3 (t) 

- |a(x)B 2 (t)] 008 3 « X T| - [|b 5 (t) - Ja 2 (t)B(t)] 

sin 3 “itj + c|pa 2 (t) + |pb 2 (t) + Je 3 ] oob» r\ 

+ Jp 3 oos3« tj + [|pa 2 (t) - |pb 2 (t)] oos(u + 2 m1 )tj 

+ [|PA 2 (T) - JpB 2 (t)] oob (u - 2**^ + |pA(x)B(t) 

sinf® + 2^)1) + |pa(x)B(x) sin^w^-a )tj + |p 2 A(x) 

oob ( u^+28 )t| + ^P 2 A(x) oob( u j-2o> )n + J-pSfT) 

sln( “j+28 )tj + Jp^x) sinC^ - 2 W )t). (J.21) 

Pram the above relation the coefficients a^ and b^ ( Introduced 
in eqns. (3*16) and (3.17)) of cob «jtj and slnujTl terms can 
be determined for various orders of subharmonio oscillations . 

(1) Subharmonio of order 2 : 

A pth harmonic component has a period 1/p of the period 
of the fundamental component. When p = where q Is an 
Integer, the order of harmonic is l/q.. Such a fractional 
order of harmonic is called a subharmonic of order q. 

Consider a subharmonio of order 2, the period of which 
is twice the period of the fundamental component. 
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Here q ■ 2, that 18 

a #/l “W /2« 

Substituting for la the eqn. (3*21) and collecting the 
coefficients of oose/2 rj and sine /2 t) terns, we have 

a l " + + & 2jl 

= |ar 2 + |p 2 a 

* Jb 5 + |a 2 b + |p 2 a 

* JhR 2 + ifrh, (3.22) 

where 

H 2 * A 2 + B 2 . 

Now the eqn. (3.18) becomes 
A(Je 2 + |iP 2 + H cos m , + y) - HB slnu^ = 0 

B(|fi 2 + fp 2 + N cos u> x + y) + HA sin ttl - 0, (3.23) 

multiplying the first equation by B and the second by A and 
subtracting 

H hi sin ^ = 0, 

which Implies 

E 2 ■ 0 for N ^ 0 f 

that is the amplitude of subbarmonio oscillation of order 2 
Is aero as long as damping is present. Therefore subharmonic 
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oscillation of order 2 can aaAnot occur In this oase. On 
the other hand for IT = 0, that Is for p *= 1.0 (zero damping) , 
we have 

H 2 (|b 2 + y + §P 2 ) - 0, 

from which 

R 2 - 0 or Jfi 2 + y + j£p 2 = 0. 

Ruling out the first solution (which leads to trivial 
solution) the require condition Is 

+ |p 2 + y - 0. (3.%) 

She eqn. (3« 24) gives the relationship hetveen the 
amplitude of subharmonio of order 2 and the excitation 
frequency a . 

Again from the eqn. (3.24) 

R 2 «-|[^ 2 + r ] ( 3 . 25 ) 

(ii) Subharmonio of order 5 t 
For this case Wj. = w /3. 

Substituting a ■ “/3 in (5*21) and collecting the 
coefficients of cos e/3 rj and sin »/3 i) v we obtain 

a l « |a 5 + |AB 2 + Jp 2 A + |pA 2 - |pB 2 
■ Jar 2 + |AP 2 + |p(A 2 - B 2 ) 
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1>1 m + ^BP 2 - felB, 

with these values of a;L and b r the eqp. (3.18) becomes 
+ |? 2 + B ooe “ a + r] + |p(A 2 - B 2 ) - KB slnSj^ . 

B[|a 2 + $? 2 + B oos » x + Y ] _ | PAB + Hi sin ^ . 0 . 

let S * j^EL 2 + jjp 2 +• K coso^ + y 

03iea the eqns. (3.26) and (3.27) become 


0 

(3.26) 

(3.27) 


AS — KB ala s . | p(i 2 - B 2 ) 

BS + KA Bln « x « | PABs 
Squaring and adding 

S 2 + B 2 sin 2 = |p pV. 
Substituting for S and simplifying, we hare 


^1 + E§* 2 + i( ff oob + y)] + 3 Y p 2 + y 2 

+ H(K + ooBu^Zf + 3P 2 ) ) ] = c, (3.29) 

where 

% - a 2 

P * P/(2 oos o + o)« 


for H = 0, that la for dissipation free system 


fe*l + <§ p2 + Jr) % + (|p 2 + y) 2 - 


0 . 


(3.30) 
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(ill) Subharmonlo of order 4 : 

\ = “/q. * »/4. 

Then from the eqn. (3.21), we have 

a-j^ ■ |a 2 + |ab 2 + |pa 2 = |ah 2 + |ap 2 
t»i = |b 3 + - |br 2 + |hp 2 . 

Combining the values of a x and b x with eqn. (3.18) 

^AH 2 + |ap 2 + H(A oos - B slnu x ) + y A = 0 

|bR 2 + |p 2 A + H(A sin “j. + B oos u^) + y B « 0, 

that is 

A[|a 2 + |p 2 + H oos + y] - KB sin w x = 0 (3.31) 

B[|a 2 + |p 2 + N oos « x + y] + NA sin w x = 0. (3.32) 

Multiplying (3.31) by B and (3.32) by A and subtracting. 

H 2 !! sin B 0» 

which implies that the amplitude of 4th order subharmonlo 
oscillation is zero whenever there is damping in the system. 
This shows the subharmonic oscillation of order 4 cannot 
occur in this case. However, with zero damping, the eqns . 
(3.31) and (3*32) reduce to 
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fa 2 + fp 2 + t - o, (3.35) 

where uy = 2 cos «/4 + a« 

(±v) Subharmonlc of order 5 i 

With “j = u /5 f following the steps given above* we have* 

|a 2 +£ 2 + Y«o f (3.34) 

where 

|iy ■ 2 cos e/5 + a. 

In a similar way the subharmonio of any order can be 
Investigated. 

Pig. 3*1 and 3.2 show the amplitude-frequency characterir* 
stic curves for different values of input amplitude P for 
subharmonic oscillation of order 2 and 5 respectively (without 
damping). The analysis reported in this chapter is valid only 
for small, detuning (y of order unity)* that is the frequency of 
subharmonio oscillation of order q is nearly equal to the 
natural frequency of the system. Sow a =-1.0, corresponds to 
a natural frequency of the linear oscillatory system* 

« o = a/3. Thus the response curves In the Pigs. 3.1 and 3*2 
can be usefully used only for u in the neighbourhood of 2nc/3 
(Pig* 3*1) and 5w/3 (Pig* 3*2). A more complete analysis for 
« * quite different from 2ic/3 or 5u/3 oan be oarried out 

using large detuning theory* which, however* is yet to be 
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developed for discrete time systems. 

Pig. 3.3 shove the variation of 3rd order suhharmonio 
amplitude with the input frequency both with and without 
damping. In this case* for each value of frequency to there 
are tvo values for the amplitude B^. As in the case of 
continuous time systems it is perhaps possible to prove that 
one of the amplitudes is stable and the other is unstable. 
However, the stability analysis results in Hill type equations, 
whose analysis has not yet been developed for discrete time 
systems. 

Next we will consider the possible sup erharmonio oscilla- 
tions in this particular case. 

(v) Sup erharmonio of order 2 t 

Here ■ 2 «• • 

Substituting this relation in (3.21) and collecting the 
coefficients of oos 2 * t) and sin 2« tj terms, we obtain 

a x - a 2 + § P 2 ) 
bj - B(jf a 2 + § p 2 ). 

With these values for a^ and b^ eqn. (3.18) takes the following 
form s 

A(|r 2 + |p 2 +y + Noos« 1 ) - NB sin » x = 0 
B(fe 2 + |P 2 + T + N oos » 1 ) + HA sin - 0, 
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from which., as before, there is no superharmonio os dilation 
of order 2 when the damping is nonzero. Whereas second 
order superharmonio solution is possible if IT ■ 0, that is 
P =1.0, Under this condition 

|a 2 + $P 2 + T - 0, (3.35) 

where jiy ■ 2 cos 2 w + a. 

The variation of amplitude with input frequency given by 
eqn. (3.25) is plotted in Fig. 3.4. 

(vi) Superharmonio of order 3 t 

In a similar manner as above 

*1 - A(|r 2 + |p 2 ) + Js 3 

- B(|fi 2 + |p 2 ). 

Then from eqn* (3*18) 

A(|fi 2 + JjP 2 + IT oos + y) + jp 3 - BB sin - 0 (3*36) 

B(^El 2 + |p 2 + V oos + y) + Sk sin* x * 0. (3*37) 

Further 

AS - NB sin | P 3 

BS + BA sin ^>0, 

3 ■ |a 2 + |p 2 + S oos ® 1 + y# 


where 




u> 


FlO-3-4 AMPLITUDE CHARACTERISTIC FOR SUPER H AR MONIC 
OF ORDER 2 /. 
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Squaring and adding the above equations, we have 
H 2 (S 2 + H 2 sin 2 » 3 .) - • 

Then substituting for S, we obtain 

+ <|p 2 + \ ff oos + |y)B| + [^P 4 + 3ST 2 + y 2 + 3?^ 

+ H oos « x (2y + 3P 2 )]B 1 - fg - 0, (3.38) 

where 

. *1 “ * 2 

and 

■ 3« • 

The variation of the amplitude of superharmonlo of order . i 
3 Is given in the equation (3.38). 

For dissipation free system, that is for N ** 0, 

(ft - 1.0) 

fgsj + ($e 2 + fr) + ($* 4 + 3*** + y 2 ) \ - §j - o (3.39) 

The results obtained for various sub/super harmonic 
solutions using the discrete tine perturbatlonal technique 
sure verified using 'Harmonic Balance* method in the following 
section. 

3.5 Method of Harmonic Balance s 


In continuous time system this is a method of wide 
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utility for obtaining periodic solution of a nonlinear 
differential equation. The study of super and subharmonie 
oscillations has also been carried out by various authors 
including Eayashi [4]» Stoker [5] and Eh [7] adopting this 
technique. The very sane technique is used in the present 
study to investigate the possible sub/super harmonic solutions 
in a strongly excited weakly nonlinear second order difference 
equations with polynomial nonlinearity. The model under study 
is same as that given in oqn, (3.1), that is 

x(k+l) + ax(k) 4 - px(k-l) + iif(x(k), x(k-l)) = E cos u k 

(3.40) 

Now the procedure for obtaining the subharmonie solu- 
tions is given below s 

Here the forcing frequency is assumed to be p* instead 
of o in the eqn. (3.40) for convenience. It will be noticed 
that, since the frequency of the external force is pu , the 
subharmonie components of order p has a frequency u and it 
may be expressed as a linear combination of cosm k and slmu k. 
Now the system eqn. (3.40) can be rewritten as s 

x(k+l) + ox(k) + px(k-l) + pf(x(k), x(k-l))= F cos p» k. 

(3.41) 

The subharmonie solutions can be developed by assuming 
the solution x(k) to the eqn. (5.41) of the form 
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x(k) « A ooa tt k + B sinu k + P oos p « k f (3*42) 

where 

A y B are amplitude functions of subha rmonic oscillation 
of order p and P is the amplitude of the fundamental compo- 
nent. as in continuous time systems, the amplitude P may be 
approximated by 

p _ F 

x “ cos pu + oj * 

Then from eqn. (3.42) 

x(k+l) = (A oos u + B sinjs ) cos o k - (A sinu - B cose ) 

sinu k + P cos p » . cos p u k - P sin p u • sin pu k. 

x(k»l) a (A cos u - B sinu ) cos u k + (A sin u + B cos « ) 

sin a k + P cos p » • oos p u k + P sin pu # sin p « k. 

Substituting for x(k+l), x(k) and x(k-l) in the eqn* (3*41), 
we have 

(A oostt + B sin « ) oos u k - (A sinu - B oosu ) sinu k 
+ P cos p u • oos p ojk - P sin pu • sin pu k 

+ oA oos «) k + aB sin u k + aP oos po k 

+ 0(A cos tt - B sinu ) cos u k + 0(A sinu + B oosu ) 

sin u k + 0 P oos p u • oos p « k + 0P .. . 

sin p « . sin p tt k + pa^ oos o k + pbp sin u k 

+ p'ag oos po k ■ F oos pu k, 
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where 

a p , D p and a 2 are the coefficients of cos « k, 
sin u k and cos pu k terms In Fourier series expansion of 
the nonlinear funotion In the eqn. (3*41) . 

Then equating the terms involving cos“k, sin to k 
separately to sero, we obtain 

cos « k t 

A cos u + B sin « + aL + p (A cosu - p sin » ) + jia^ = 0 

A(p+1) cosu + B(l-p) slnu + aA + |i a p ■ 0. (3.43) 

sin to k : 

-A sin tt + B cos m + aB + p(A slnu + B cosu ) + |ibp * 0 

B(p+1) cos u + A(p— 1) sin u + aB + [i bp = 0* (3.44) 

The eqns. (3.43) and (3.44) give the variation of sub- 
harmonic amplitude with the input frequency. The method of 
obtaining the amplitude of various orders of subharmonic 
oscillations are illustrated by considering a particular 
nonlinearity as indicated In the following example. 

3.6 Example : 

Consider the following nonlinear difference equation 

x(k+l) + ax(k) + px(k-l) + |xx^(k) = F cos p u k. (3*43) 

Eqn. (3.45) is same as the eqn. (3.41) except for the nonlinear 
funotion, that is given by 
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ttO-x’OO. 

Then the interest is to find explicit expressions for a^ f 
bp and a 2 « 

Jw » a oostt k + B sins k + P cos pw k) 3 
*= A 3 oob 3 s k + 3 A 2 ! cos 2 w k. sins k 

+ 3 AB 2 008 wk * sin » k 
+ B 3 sin 3 « k + 3P oos puk (A 2 cos 2 u k 

2 2 

+ 2 AB ooStf k. sin a k + B sin u k) 

+ 3P ? oos 2 p «k (A oos« k + B sins k) + P 3 oos 3 p «k» 
Expanding and simplifying 

* 5 (k) - (f A 3 .+ |iB 2 + |«lP 2 ) oos u k + (| A^ + |b 3 + \ BP 2 ) 
sin u k + ( J A 3 - | AB 2 ) cos 3u k + (J- A% - J B 3 ) 
sin 3* k + (| A 2 p + \ bS + |p 3 ) oos p w k + jj$ 3 oos 3 Jfc k 
+ (J A^P - Jb 2 ?) oos(p u + 2b )k + (|a% - |iB^P) 
oos(p« - 2u )k + |&BP sin(2u +■ pu ) k + |iBP 
sin(2s - p« )k + ^AP 2 cos(u + 2pi )k + |&P 2 
oos ( co - 2pe )k + |fiP 2 sinCu 4- 2po» )k + |fiP 2 


sin (» - 2p s)k. 


(3*46) 
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Now using the eqn. (3.46) the various orders of sub- 
harmonic oscillations are studied as follows. 

(i) Subharmonic of order 2 : 

Here p ■ 2. t 

Then from the eqn. (3.46) 

a,, « | A 3 + J AB 2 + | AP 2 
= | AH 2 + \ AP 2 

b 2 s + I® 5 + \ BP 2 

= |hR 2 + | BP 2 , 

where 

R 2 - A 2 ♦ B 2 . 

Substituting a 2 and b 2 in oqns. (3.43) and (3.44), we 
obtain 

A(p+1) oob w - B(p -1) sine + oA + ji 4 (Jr 2 + |p 2 )= 0 (3.47) 

B(p+1) oos o + A(p-l) sins + op + (iB(|& 2 + |p 2 )= 0, (3.48) 

that is, 

A[ (p+1) cosu>+ a + |i(jj H 2 + |iP 2 )] - B(p-l) sin co ■ 0 

B[ (P+1) oob « + a + ji(J- H 2 + !jp 2 )] + A (P-1) sin u « 0. 

Multiplying the first equation by B and the second by A and 
subtracting. 
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0-1) R 2 slnu = 0, 

that Is R = 0 for (1-p) 0, which implies the amplitude of 

subharmonic oscillation is zero as long as the damping is 
present In the system. However with zero Aar *p~i (0 = 1.0) 
we hdLve, R 4 0 and aubharmonic of order 2 can occur, that is 

2 cos « + a + \i (Jr 2 + J P 2 ) s* 0* 

Pram which 

| E 2 + |f Z + y - 0, (3.49) 

where 



and 

PY ** 2 cos » + a.. 

It is interesting to note that the eqn. (3.49) and (3.24) 
are identical. In (3.24) the frequency = w /2 is the 
frequenter of subharmoziio oscillation of order 2 j whereas in 
(3*49), the frequency to is the frequency of subhormonio oscilla- 
tion of order 2 (the input frequencies are o> and 2 to respecti- 
vely). 

(ii) Subharmonic of order 3 : 

Here p ■ 3. 


Then as before 
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a 5 « | A 3 + | AB 2 + | AP 2 + - Jb 2 ? 

JiR 2 + \ AP 2 + |p (A 2 - B 2 ) 
bj - | BR 2 + §BP 2 - \ ^BP. 

Combining these values with eqns. (3*43) and (3.44), ve have 

A(p+1) oos u . B(p-l) sin » + aA + |iA(|fi 2 + §P 2 ) 

+ | »iP(A 2 - B 2 ) * 0 (3.50) 

B(p+1) cosw + A(p-l) sins .+ oB + jiB(Jft 2 + |p 2 ) 

- \ pAHP * 0. (3.51) 

Xhe eqn. (3.50) can be rewritten as 

A ^(P+1) ooS u - 2 oobm + 2 OOBtt +_g + | a 2 + 

- Bina + | P(A 2 -B 2 ) = 0 
A[I*L22£JL-Jt^l + Ifizil oos u + \ B 2 + |p 2 ] 

- SUbil sin « + | P(A 2 - B 2 ) = 0, 
that is 

A(|h 2 + j^P 2 + fl cos CO + Y) + | P(A 2 -B 2 ) - BBT sinu ■ 0 

(3.52) 
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and similarly from eqp« (3*31) » we have 

B(| E 2 + |p 2 + N cos « + y) - | ABP + M sin w » 0, (3.53) 

where 

HY « (2 oos a + a), 
and N = . 

The eqns. (3.52) and (3.33) are identical to eqns. (3.26) and 
(3.27) respectively. Following the same procedure as in eqns. 
(3,26) and (3.27) we have the following equation identical to 
the eqn» (3.29), 

Is®! + p2 + 1 008 w+ Y)] h + [f p4 + + y 2 

+ Hf(H + cos to (2y + 3P 2 ))] » 0.0 

where 

^ - E 2 

and P «* F/((p>*tL) oos 3 8 - 1 - a]. 

(lii) Subharmonic of order 3 : 

Proceeding as before with P = 5, we have 

“5 -? Aa2 + l “ 2 
t> 5 - | BR 2 + \ HP 2 , 

that is with p 1.0 (with damping), we have 
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A((p+1) oos u + a + B 2 + \ E 2 )] - B(P-l) sin M ■ 0 

B[(Btl) oos at + a + |»(J E 2 + 3[ E 2 )] + A(p-l) sin « « 0, 
from which. , as before for p ■» 1.0 

J a 2 + \ P 2 + T - 0. (3.54) 

where 

E *■ F/(2 oos 5 « + a). 

Here again the eqn. (3.54) is identical to equation 
obtained by multiple scale perturbation technique. 

As in the case of multiple time perturbation technique , 
the identification of super harmonic components by harmonic 
balance method is not straight forward, but the system 
equation is to be changed slightly. She analysis is carried 
out as follows. 

Let the system equation be of the form 
x(k+l) + <xx(k) + px(k-l) + |if (x(k) # x(k-l))« E oos « k f (3.55) 

where 1 <* 1 is the input frequency. It ie of Interest to 
investigate superharmonic solution of order p (p > 1) Md so 
the solution is assumed to be in the following form s 

x(k) a A oos pule + B sin pu k + E oos to k, (3.56) 

where 

A and B are the amplitude of superharmonio solution of 


order V. 
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Then substituting the eqn. (3*56) in (3.55) and simpli- 
fying, we obtain 

A(p+1) oos p« - B(p-l) sin pu + <xA + \x a p ■ 0 (3.57) 

B(p+1) oos p» + A(p-l) sin pm + aB + p b p » 0 • (3.58) 

As in the aase of subharmonic solution ,P, the amplitude 
of the fundamental component is approximated by 

1 " T ^lT ' ooi»+" ' a • 

where 

ap and bp are the coefficients of oos pm k and sin pu k 
terms in the Fourier series expansion of the nonlinear fun- 
ction f(. f .) in (3.55) • 

let f(. ) » x^(k). ■ Vith this nonlinearity the expression 
given (3.46) ftcan be used to determine the terms a^ and bp 
simply by replacing* a and pu * « • Vith this 

construction the following superharmonic solutions ore 
studied. 

(iv) Superharmonio of order 2 : 

Superharmonic os oillation of order 2 is obtained by 
assuming p ■ 2 in the eqn. (3*46) after incorporating the 
above replacement. That is, 

a 2 - | AH 2 + | AP 2 
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b 2 * | BE 2 + | BP 2 . 

Combining these values with eqns. (3.57) and (3.58) 

A(p+1) cos poi - B(p-l) sin pe + oA + |iA(J H 2 + Jjp 2 ) - q 

B(P+1) oos p» + A(p-i) sin p M + <xB + pB(| H 2 + |p 2 ) = 0 , 
from which, as in the above cases, we have 


0 ■ 1.0 

i ® 2 + i p2 + r - 0. 

where 

HY * 2 cos 2e + a 

and 

1 

® ■ ®/(2 oo Su + a). 


(3.59) 


A® expected the e<jn. (3.59) ie identical to the ep 
(3.35). 

(▼) Superhaxmonlo of order 3 : 

Vlth p ■ 3» we have 
a 3 = | Afi 2 + | ap 2 + ^ P 3 

b 3 - \ BE 2 + \ BP 2 . 

Combining these values with eqns. (3.57) and (3.58) 
A(p+1) oos Pu “ B(P-l) sin pot + aA + |iA(^fi 2 + jjp 2 ) 

+ J » i P 3 - 0 


(3.60) 
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B(p+1) oos p« + A (P-1) sin bd + aB + iiB(|-R 2 + |p 2 ) = 0 V 

( 3 . 61 ) 

adding and subtracting 2A cos in the eqn. (3.60) and 
2B oos pa in the eqn. (3.61), we have 

A(|a 2 + fp 2 + H oos pet + y) + - BB sin pu « 0 (3.62) 

B(Jr 2 + |p 2 + N cos pa, + y ) + KA sin b„ * 0. (3.63) 

The e^ps.(3.62) and (3.63) oan be oompared with eqns. (3.36) 
and (3.37). 

Then from eqns. (3.62) and (3.63) , we obtain 

fg + (| P 2 + | H oos pw + |y) R 2 + [| B 4 + H 2 + y 2 + 3?^ 

+ B oos po, (2 y + 3P 2 )]B 1 - - 0.0 (3.64) 

where 

R^ * R 2 

and hy * 2 oos 3w + a. 

In a similar way the other superhamnonio oscillations 
can be investigated. 

3.7 Bis cuss ion of Results : 

The following observations oan be made from the results 
obtained in the previous sections for cubic nonlinearity. 

(i) An observation of the results obtained for various orders 
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of sub /superharmonic oscillations by (a) Discrete tine per- 
turbations!. technique and (b) Harmonic balance method, shove 
excellent agreement between the two methods* 

(11) The conditions deduced for the existence of sub/super- 
harmonic are shown In Figs* 3*1 to 3*4* Figs. 3.1 and 3*2 
show the variation of the amplitude of subharmonio solution 
of order 2 and 5 with input frequency. It must be pointed out 
that the existence of sub/superharmonic oscillations has been 
established only for the values of detuning parameter y of 
0(1), the above curves con not effectively be utilised for the 
following reason. As mentioned earlier, the natural frequency 
« 0 of the linear system is w/3, and the value of r is of order 
unity only when « ■ 2ic/3 for order 2 and u « 5n/3 for order 

5 subharmonio oscillations* She plots in Figs* 3*1 and 3.2 
show that value of the detuning factor y works out to be very 
high hence with the present approach, small detuning theory 
cannot be used for detecting order 2 and order 5 subharmonic 
oscillations* The same argument is applicable for the super- 
harmonic oscillation of order 2 with reference to the Fig* 3*4* 

(Hi) For subharmonio oscillation of order 3, the variation of 
amplitude with input frequency is shown in the Fig. 3*3 with 
zero and weak damping in the system. Here again the natural 
frequency « 0 of the linear system is «/3# For this sub- 


harmonio oscillation the value of the detuning parameter y 
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Is expected to be of order unity only when the Input frequency 
u ■ w. For F * 1.0 as veil as for F = 3.0 , with damping, 
large detuning theory is necessary for effective utilisation 
of the response curves. However, with zero damping it is 
dear from the figure, small detuning theory analysis is 
possible and hence a subharmonio oscillation of order 3 con be 
sustained. As in the case of continuous time systems £4,117] 

the sub/superhaxmonio oscillations are highly dependent on 

■ 

Initial conditions and furthermore sinoe the equilibrium 
points in the A(x) - B(x) space are expected to be centers/ 
saddles (zero damping), henoe the time response of the given 
system shoving the subharmonio os oil. Tat ion of order 3 is all 
the more difficult to obtain. 

(iv) It is also to be pointed out that for the other values 
of F, namely F=2.0, 4.0, etc. (not shown in the Fig, 3.3)# the 
structure of the curves are such that the impressed frequency, 
for admissible amplitude is further away from * than for 

F * 3.0. 

(v) Another interesting observation is, that, for a given 
input amplitude (F = 1.0), the curve shown in the Fig. 3-3 is 
not dosed as indicated by points A and B. It is expected 
that the two branches would nest, but this would entail finer 
increments in “ than what has been considered. In any case 
the dosed curves would still be not in neighbourhood of * for 
effective use. 
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5*8 . Conclusion ; 

The problem of sub/sup erharmonlc oscillations In 
nonlinear discrete time systems has been studied In this 
chapter by means of discrete multiple soale perturbatioual 
technique as well as by harmonic balance method and conditions 
for existence of various orders of sub/superhaxmonio oscilla- 
tions have been established. An example was presented with 
cubic nonlinearity and analysed In detail to demonstrate the 
use of the derived conditions. Proper Initial conditions 
have, however, be selected to actually see the sub/super- 
harmonio oscillations in the time response. 



CHAPTER 4 


PERIODIC OSCILLATIONS IN DIGITAL FILTERS 
4-1 Introduction i 

In the earlier chapters a class of second order nonlinear 
difference equations with and without external input were con- 
sidered. Polynomial kind of nonlinearity was assumed and it 
was demonstrated that the multiple tine perturbation technique 
can be applied to study the system behaviour. In this chapter, 
a second order difference equation with a nonlinearity other 
than the polynomial type is considered. As mentioned in the 
introductory chapter the description of operation of digital 
filters due to finite wordlength registers typically involves 
nonlinearities other than those of polynomial type. The 
various possible nonlinear effects that can occur in digital 
filters due to the finite wordlength available for the signal 
storage are briefly explained in the following section. 

4.2 Nonlinear Phenomena in Digital Filters: 

The wordlength of the signals in a digital system, will 
' in general Increase after arithmetical operations (addition 
and multiplication) . When storing signals after arithmetical 
operations the wordlength has to be reduced which is equivalent 
to nonlinear operations being performed on the signals. Two 
types of nonlinear phenomena are observed in digital filters : 
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(1) quantization and (2) overflow. Quantization la a word- 
length reduotlon In whloh the least significant bits are 
affected £79]. It will lead to an error that la proportional 
to the value corresponding with the least significant bit of 
the quantized signal. As suoh f the error will be quite small, 
but an error of this kind will be Introduced during every 
sample period and the combined effect of all these errors may 
cause a severe distortion of the filter behaviour. Overflow 
occurs in digital filters, if after arithmetical operations, 
a word results that lies outside the range of the speoified 
wordlength of the register. In this oase the nost significant 
bits of the signals must be affected before it can be stored 
in the register and this nonlinear phenomenon introduces large 
errors . 

Both the quantization and overflow nonlinearities will 
sometimes generate sustained oscillations known as limit cycle 
oscillations in digital filters. Limit cycle oscillations 
due to quantization have relatively small amplitudes and they, 
oannot be eliminated so easily. They have been studied In 
considerable detail using fixed point arithmetic [43, 79, 

93 , HI] • On the other hand the overflow nonlinear operation 
introduoes a very large error and large amplitude limit cycle 
oscillations are then possible. Methods have been reported 
for analysis of these effects [76, 79, 64] and to suppress 
them to certain extent [96, 97] • A- 8 mentioned earlier, the 
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effects due to overflow nonlinearitios aro Quite severe end 
their effects are analysed In detail In the present study. 

The influence of the quantisation effects will assumed to be 
small and will not be considered in this analysis. 

She study of overflow behaviour manifesting nonlinear 
oscillations because of adder overflow, and oc curing in a 
second order section of a digital filter was begun in [76] 
and [80].. The system considered is Shown in Fig.. 4.1 and 
the various possible overflow nonlinearities are given in 
Fig. 4.2. It was shown in [76] that for a certain type of 
overflow nonlinearity; for example twos complement shown in 
Fig. 4..2(c), limit cycle oscillations can be sustained, for 
appropriately chosen initial conditions, under zero input 
situation. It was shown there that such oscillations can not 
exist in the zero input response of a seoond order digital 
filter employing saturation arithmetic depicted in Fig. 4.2(a). 
However f such limit cycle oscillations are possible with satu- 
ration nonlinearity if one takes, the quantization effects also 
into account. Willson [84] among others has analysed a seoond 
order digital filter under zero input condition considering 
quantization effects. This analysis is most general and can 
be easily extended to wide varieties of overflow nonlinearities. 
Recently Mitra [92] made on observation that amongst digital 
filters employing the saturation arithmetic a fundamental diff- 
erence exists between seoond order and higher order sections. 
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sanely the latter may sustain overflow oscillations. In all 
the above analysis it was assumed that the response of the 
linear model Is asymptotically stable, that is the filter 
coefficients a and b have values determined by points that 
He within the open triangular rlglon shown In Fig. 4.5. 

An excellent review on effects of finite wordlength Is given 

by Liu [40] and the limit cycle problem Is given In [79, 98]. 

\ 

Xhe available results on limit cycles concentrate 
on, conditions for the existence of limit cycles, estimate of 
limit cycle amplitude and frequency. A bulk of work has 
concentrated an second order digital filter with coefficient 
values Inside the stability triangle. There has however been 
no work on limit cycle os dilations in second order digital 
filters for the coefficient values outside the stability 
triangle ABC given in Fig. 4.5. It is therefore of Interest 
to study 'the limit cycle oscillations outside the stability 
triangle under foroe free situation. In this Chapter, the study 
of overflow oscillations is extended for the parameter values 
outside the stability triangle with saturation overflow non- 
linearity. 

She following are the main contributions in this study : 

(1) The existence of limit cycles with different periods is 
demonstrated for parameter values outside the stability 
triangle. Regions delineating the parameter space 
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have been obtained whore a particular kind of Unit 
cycle os dilation with a known poriod is obtained. 

(li) To. investigate a region in the parameter plane in which 
the output sequence decays monotonically for a set of 
values in the initial condition plane. 

4.3 Limit Qy ole Os dilation Outside the Stability Triangle t 

It is evident that for parameter values outside the 
stability triangle* the linear response of the filter is un- 
stable, which in turn results in nonlinear oscillations due to 
saturation nonlinearity. This nonlinear oscillation ig however 
not always present for some values of £a,b] outside the stability 
triangle as demonstrated later In the analysis. In this section, 
the study is concentrated on identification of various regions 
outside the stability triangle, in the parameter plane where 
Unit cycle oscillations of different periods are possible. 

The system under Investigation is depicted in Fig. 4.4 and its 
mathematical description is given as 

z(kf-2) * f[ax(k+l) + bx(k)] (4.1) 

where a and b are the filter coefficients or paraaeters and 
f ( •) is the nonlinear function representing a saturation non- 
linearity described by 
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1.0 

1.0 (4*2) 

1*0 

with the schematic representation of this saturation non-!* 
linearity shown in Fig. 4.2(a). 

4.4 Analysis t 

It is well known that within the stability triangle in 
the a-b parameter space the free system described by 

x(k+2) = ax(k+l) + bx(k) 

has an asymptotically decaying form for all initial conditions. 
Of interest is the nature of the solution for parameter loca- 
tions outside the stability triangle and In particular the 
identification of regions where a limit cycle oscillation of 
given period L (denoted by P^) is obtained. Three methods 
have been proposed for the above purpose, namely 

(a) Baploylng actual saturation nonlinearity (Fig. 4. 2a) 

(b) Using an approximate polynomial for saturation nonlinearity 
and employing harmonic balanoe method. 

(c) Using an approximate polynomial for saturation nonlinearity 
and applying disarete multiple scale perturbatisnal 
technique [100]. 
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4.4.1 Aaploying actual saturation nonlinearity i 

The technique proposed by Ebert et al [76] for location 
of regions of P^ and ?2 limit cycle oscillations vlth twos 
complement arithmetic (Fig. 4 • 2d) has been extended for the 
present study , that is for investigating regions outside the 
stability triangle in which P-^ limit cycle oscillations are 
possible employing saturation overflow nonlinearity. Thus 
considering the second order digital filter given in eqn. 
(4.1) and the nonlinearity given in eqn. (4.2) the following 
P^ oscillations are discussed* 

P^ Limit cycle oscillations 

In this case the output of the filter x(fc) is a constant 

[76] 

Let x(k) - A 
Substituting In eqn. (4.1) 

A » f[A(a + b)] 

Then using the definition of the saturation nonlinearity 
described in eqn. (4*2). 

If [A(a+b)] >■ 1.0 

A = 1.0 (4*3) 


If 


[A(a4-b)] < - 1.0 
A = -1.0 


(4.4) 



and if I A(a+b)l <. l.o 
A * A (a+b) 

i.e. a+b ■ 1.0 (4*5) 

the above expression in (4*5) represents the linear operation 
of the filter. 

Pron eqns. (4*3) and (4-4 ) $ we obtain 

a + b > 1.0 (4.6) 

The region in the a-b plane where a constant output of 
magnitude unity exists is given by the inequattfy (4.6). 

P 2 Unit cycle oscillation : 

In general the period two periodic oscillation is 
expressed as 

x(k) * A + B cos nk (4.7) 

where, 

with A 4 0 and B 0, the oscillation is about a con- 
stant d-c value > whereas with A » 0 and B ^ 0 the oscillation 
has no d-c bias. Then from eqn. (4.7) 

x(0) * A + B 
x(l) *sA “ B 

and this sequence (A + B, A - B) repeats, for all even and odd 
values of k respectively. 
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Substitution of the values x(0) and x(l) In ean, (4*1) 

leads 


A + B - f [a(A-3) +.b(A+B)] 
A - B - ffa(AtB) + b(A-B) ] 

fron the above 

if |>(A-B) + b(A+B)] > 1.0 

A + B - 1.0 

if (>(A+B) + b(A-B) ] > 1.0 

A - B = 1.0 

if [a(A+B) + b(A-B)] < -1.0 
A - B * -1,0 

if [a (A-B) .+ b(A+B)] < -1,0 
A + B - -1.0 

whereas 

if | a(A+B) + b(A-B)| < 1.0 

A-B « a(A+B) + b(A-B) 

and 

if | a (A-B) + b (A+B)J < 1.0 
A + B * a(A-B) + b(A+B) 

Finally 

if [a(A+B) + b(A-B)] « 0 
A - B - 0 

if [a (A-B) + b(A+B)] » 0 


(4.8) 

(4.9) 

(4.10) 
(4,11.) 

(4.12) 

(4.13) 

(4.14) 


A + B 5 0 




(4,19) 
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Vow the required region in the a-b plane can be obtained 
by considering the above equations as follows : 

Solving eqns. (4*8) and (4*9) we get, 

A * 1.0 
B - 0.0 
with a+b > 1.0 

and si m il ar ly considering eqns. (4*10) and (4.11) 

A = -1.0 
B * 0.0 
with a+b= > 1.0 

which in fact shows that P^ Unit cycle oscillation of unit 
magnitude is possible In the region a+b > 1.0 which is the 
same result obtained In eqn. (4.6), 

Further considering eqns. (4.8) and (4.10), the following 
relations are obtained 

A » 0.0 
B - 1.0 
and b-a > 1.0 

and In a similar way with equations (4.9) and (4*11) we 
obtain y 

A = 0.0 


B * -1.0 
and b-a > 1.0 
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which implies that limit oyole oscillations of period two 
without d.o. component is possible in the region 

b - a > 1.0 . (4.16) 

Further by combining eqns. (4.8) and (4.14) we get 

A * 0,5 
B = 0.5 
a = 0,0 
and b > 1.0 

Similarly the following conditions are obtained by 
considering eqns. (4.9) and (4.15) 

A - 0.5 
B —0.5 
a * 0,0 
b > 1.0 

Also from eqns. (4.10) and (4.16) 

A —0.5 
B * 0.5 
a = 0.0 
b > 1.0 

further from eqns. (4,11) and (4.14) 

A —0.5 
B —0.5 
a ■ o.o 
b > 1.0 
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which, shows limit cycle oscillation of period two with d.o, 
bias, of the form x(k) * + 0.6 + 0.5 ooswk is possible only 
on the line b > 1.0 and a * 0,0. 

Pj limit cycle oscillation t 

In this case the general solution is assumed to be of 
the form 

x(k) = A + B oos + c sin 
then 

x(o) ■ A + B 

x(l) * A - B/2 + o 

x( 2) - A - B/2 - ^ o 

and this sequence repeats for all k > 3. 

Then the substitution of the above relations in eqn. 
(4.1) leads the following equations, 

x(2) ■ A - | - ^|o = f[a(A - | ^ o) + b(A + B)] 

x(3) - A + B ■ f[a(A - B/2 - 4|o) + b(A - B/2 + o )] 

x(4) - A - f + *|e « f[a(A+B) + b(A - B/2 - %>)] 

if 

a(A - | + o) + b(A + B) > 1,0 
A - ^ ^ c * 1*0 


(4.17) 



IS 5 

If «(A - | - Jt| p) + b(A - | + ^| o) > l.o 


A + B - 1.0 

(4.18) 

If a(A+B) + b(A “ | ^ c) > 1,0 


A - | + o « 1.0 

. (4.19) 


Solving (4*17)* (4*18) a n d (4*19) wo get 

A * 1.0 
B ■ 0.0 
0 = 0.0 

that is the solution x(k) * l.o and the corresponding 
regin la 

a + b > l.o 

which Is same as (4.6). 

Let 

a(A - f + o) + b(A + B) b o.O 

A - | o - 0.0 (4.20) 

and 

a(A+fi) + b(A - | ^ o) < •l.o 

A - ! + o - - 1.0 . 

Solving eons. (4.18), (4.20) and .(4.21) we obtain 

A = 0.0 
B * 1.0 



(4.21) 



166 


and the corresponding regin In the a - b plane Is 

a * b 
b < - 1.0 

a < - 1.0 . (4*22) ■ 

Bln. (4.22) gives the straight line, on which limit cycle 
oscillation of period 3 is possible without d.o. component. 

Further for 

I a(A - | + o) + b(A + B)| <. 1.0 

A ~ § " ° " a ( A - § + c) + b(A + B) (4*23) 

then the oo editions given In (4.18 ) » (4*21) and (4*23) lead 
to the following region There a limit cycle oscillation of 
period 3 with d.c. component can be sustained. 

I b - a| < 1.0 

a(b-a) - b> 1.0 

a + b(b-a)< -1.0 (4.24) 

limit cycle oscillation : 

The output sequence takes the form 
z(k) ■ A + B cos jS + o sin 


then 
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x(o) = A + B 
x(l) - A + 0 
x(2) * A - B 
x(3) = A - C 

and this sequence repeats for all k £ 4. 

Then, 

x(2) « A - B * £[a(A+C) + b(A+B)] 

x(3) - A - C a f[a(A-B) + b(A+0)] 

x(4) - A + B - f[a(A-C) + b(A-B)] 

x(5) a* A + 0 as f[a(A+B) + b(A-C)] 

These equations combined with the properties of the 
saturation nonlinearity, the region In which a limit cycle 
oscillation of period 4 Is possible is given by 

a + b < -1.0 

a - b > -1.0 (4.25) 

In. a similar way a region in the .a - b plane in which, 
a period 6 limit cycle oscillation can be sustained Is given 
by 

| a + b| < 1.0 

a - b(a+b) > 1.0 

a(arf-b) + b< -1.0 (4.26) 

It is Interesting to observe that the region of Pg oscillation 
Is identical to region of Pj os dilation mirrored around a =0.0. 
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that la 

f(y) « 1.16253jr - 0.285J 5 + 0.03 64 3y 5 - 0.00167Jr 7 , with 

Ul i 3.0. t;«(4*2 1 

The saturation nonlinearity and its polynomial approximation 
(dotted lines) are indicated in Pi g. 4.5. 

The approximate polynomial expression given in (4*27) can 
be used to investigate the regions in the a-b plane in which 
various limit cycle oscillations are possible. The analysis 
proceeds as follows. 

P 1 limit cycle oscillation : 

Let x(fc) * A, a constant 
Prom the system configuration shown in Pig. 4.4. 

y(k) ■ ax(k+l) + bx(k) 
y = A(a+b) • 

Again from Pig. 4 .4, 

x(k**2) = f[y(k)] 

that is A s f(y) = 1.16253y - 0.285y 3 + 0.03643y 5 - 0.00l67yj 
Substituting for y 

A * 1.16253A(a+b) - 0.285 A 5 (a+b) 3 + 0.03643A 5 (a+b) 5 

-0.00l67A 7 (a+b) 7 

with 

| A(a+b)| <3.0 

the above equation can be rewritten as ■ 
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FIG-4-5 APPROXIMATION TO SATURATION NONLINEARITY 
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1.16253 - 0.285 A 2 (a+b) 2 +0.03643A 4 (a+b) 4 .- 0.0Q167A 6 (a+b) 6 

* a+b , (4*. 28) 

Substituting the maximum bound on | A(a+b)j which Is 3.0 f 

a + b - 3.02. (4.29) 

This gives the -yi-nm bound on (a+b ) • 

In a similar way the minimum bound on (a+b) c** be 
obtained by substituting |A(a+b) J = l.o, 

a + b * 0.98 . (4.30) 

From equations (4*29) and (4.50) the region in which a limit 
cycle oscillation of period one is possible is given by 

0.98 < a + b < 3.02 # 

It is also noticed that as D booonos larger the upper bound on 
a+b also becomes greater and tho ooaplete region in the a-b 
piano for D is given by 

a + b > 0.98 (4131) 

Pg limit cycle oscillation t 

In this oase the output sequence assumes the form 

z(fc) - A + B oos itk 
y(k) ■ aac(k+l) + bx(k) 

■ A (a+b) + (b-a) B oos nk 

■ P + Q oos wk 
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where 

P • A(a+b) 

Q «B(b-a) 

Considering the approximate polynomial for . saturation 
nonlinearity 

x(k+2) * f(y) » 1.16253? - 0.285y 3 + 0.03643? 5 - 0.00167y 7 
with l y | * | P + Q cos wk | - 3*0 

that is 

A + B oos wk m 1.16233 (P+Q cos wk) - 0*285 (P + Q cos wk) 3 

+ 0.03643 (P + Q ooa wk) 5 - 0.00167 (P + Q cos wk) 7 
■ ■ 

Expanding and equating the constant term and the coefficient of 
oos wk term separately using the following trigonometric 
relations 

oos n wk = sos wk for n odd 

m 1.0 for n even 

■ ■ • " 1 

ye obtain 

P[1.16253 - 0.285 (P 2 + 3Q 2 ) + 0.03643 (P 4 + 10 P 2 Q 2 + 5 0 4 ) 

- 0.00167 (P 6 + 21 P 4 Q 2 + 35 P 2 *? 4 + 7 Q 6 ] * A 

Q[l. 16253 - 0.285 (5P 2 + Q 2 ) + 0.03643 (5P 4 + 10P 2 Q 2 + Q 4 ) 

- 0.00167 (7P 6 + 35P 4 Q 2 + ZB 2 # + Q^] ■ B m (4.32) 
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Oaae 1 

If A = 0, then P = 0 and for a ^ b 

and also |y | < 3.0 implies | Q| < 3*0.. 

Then eqn. (3.32). reduces to 

1.16253 - 0.265 Q 2 + 0.03643 Q 4 - 0.00167 Q 6 * . 

Substituting the bounds on Q 

0.98 < b - a < 3.02 . 

Then as indicated earlier the region in a-b plane (for D ■*» ) 

In which ?2 limit cycle oscillations are possible is 

b - a > 0.98 „ (4.33) 

<fe.Pg. U 

If A 0 and B * 0 

Conditions are same as for P^ limit oycle oscillation. 

Oase ill 

If A 0 and B ^ 0 

Then the equations in (4.32) can be rewritten as 

1.16253 - 0.285 (P 2 +3Q 2 ) + 0.03643 (P 4 + 10P 2 Q 2 + 5Q 4 ) - 0.00167 
(P 6 + 21P 4 Q 2 + 35P 2 Q 4 + '7Q 6 ) - (4.34) 

1.16253 - 0.285 (3P 2 +Q 2 ) + 0.03643 (5P 4 +1QP 2 qW) - 0.00167 

• ■ * 

(7P 6 +35S*Q* + ZUpV + Q S ) - ^ (4.35) 



174 


Solution to eqns* (4.34) and (4.35) la possible only if a « 0. 
Under this condition P * Q, which implies A m 3, that is 

x(k) = A(1 + cos *k) 

which is of the same, form obtained by the first, method* In 
a similar way P^,?^ and Pg limit cycle oscillations can be 
obtained. 

Xhe method of obtaining the above regions using the 
discrete multiple time perturbation technique employed in 
Chapter 2 Is mow given. 

4.4.3 Approximate polynomial and multiple scale perturbation 
method s 

A multiple time perturbation technique has been proposed 
in [100] to study a class of discrete time system with polyno- 
mial kind of nonlinearity for obtaining an approximate oloBed 
form solution. Xhe very same technique la employed here to 
Classify different regions in a-b parameter plane in which 
different Pj os dilations are possible by suitably rearranging 
the system equation vith polynomial approximation for the non- 
linearity. 

Xhe pertinent filter equation to be analysed is given by 
x(kf2) = 1.16253y - 0.2857 5 + O.Q3643y 5 - 0.0<XL67y 7 
where y ■ y(k) = a‘x(k+l) + bx(k) # 
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Xhe above equation can oonvenlnntly be written as 

x(k+2) + a 1 x(kfl) + a^x(k) * iig[x(k),x(k+l)] (4.36) 

where g is a known function of lndloated variable and is 
given by 

tig(«) * (1.16253a -i- a 1 ) x(k+l) 4- (1.16253b + a 2 ) x(k) - 0.285 
(ax(fcfl) + bac(k)] 5 + 0.03643 [ax (kfl)+bx(k)] 5 

- 0.00167[ax(k+l) + bx(k)] 7 . (4.37) 

and a^ and a 2 are constants to be selected such that the Hwiw 
response for p « 0 is periodio with known period. 

limit cycle ob cilia t ion 

From the properties of the stability triangle, for 
a^ = -a and & 2 * -p such that a + p « 1.0, the response of 
the linear model is given by 

x(k) rn A + B(-p) k . (4.58) 

Vlth the above values for a^ and a 2 eqzu (4.36) oan be 
rewritten as 

x(k+2) - ox(kfl) - px(k) ■ |ig(x(k), x(kfl) (4.39) 

where 

pg( *) - (1.16253a - «) x(ta-l) + (1.1625.3b - p) x(k) 

- 0.285[ax(k+l) + bx(k)] 3 + 0.03643 (ax(Jefl) + b^(k)] 5 

- O.OOL67 (ax(k+lj + b;x(k)] 7 # (4.40) 
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The eqn. (4.59) is a weakly nonlinear difference equation 
for which multiple time perturbation 'technique oan be applied. 
Here the dependent variable x ie a function of two Independent 
variables namely r\ and x and so the variable z and its diffe- 
rences oan be expressed In taros of tj sad x as follows. 

x(k) » x(tj,x) 

x(k+l) a x(ij+l»*0 + |ix(t],T+l) - |i x(ri,x) 

x(k+l) a x(q+2,x) + 2)1 x(t)+l, x+1) - 2)ix(t)+1,t) . (4.41) 

It is to be noted that eqn. (4.41) oan easily be obtained by 
substituting I » 4 , the forward difference operator and 
= "jo a o in eqns. (2.24) and (2.25). In this present 
study the analysis is carried out through two time scaling, 
model obtained using the forward difference operator, since 
the Interest Is to obtain conditions under steady state situa- 
tion. 

Xhe solution x(k) is also expressed In the following 
asymptotlo series as • t 

x(k) a xfa,*) a l ^.{t^x) . (4.42) 

jaO J .. 

Substituting eqns.- (4.41) and (4.42) in (4.79) and collecting 

i l 2 

the coefficients of n , |i » )i etc. terms and equating them 
separately to eero, we obtained the following set of linear 
difference equations f 
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x 0 (t|+2,t) - a x o (ti+ 1.*0 - P x 0 (n*t) -0 (4.43) 

^(yi+2,*) - a x^(t)+1»t) - p x^t^t) « a[x 0 (Tj,*+l) - x 0 (t|,O] 

+ i l s tfg j L»r *1 s Q (n+l«t) x 0 <t,^)- bf 5 - 

[a ac 0 (n+i,*p) + b XqCt ),*)] 3 
+ x 0 (u+1,t) + bx 0 (ti,T)] 5 

- [ax 0 (TH-l,T) + bx 0 (Ti f T)] 7 B (4.44) 

Solution toeqn. (4.43) Is known as generating solution and It 
Is of the form 

* 0 (t|tT) - A(t) + B(t) a (-P) 11 (4.45) 

where A(t) and B(t) are the amplitude functions and their 
explicit form can be determined by substituting (4.45) on the 
right hand side of eon. (4.44) and suppressing secular terms ^ 
Following the procedure given In Chapter 2 suppression of the 
secular terms leads to 

at AM + (1.16353a - «> a(t) + (1.16855b - BJ i(t) 

A 1 

. A 5 ( t) (a + b)5 - kfifflSZ 

A 7 (t) (a+b) 7 ■ 0 

under steady state tho above aquation takes the following form 
1.16253 7 0.285 A 2 (a+b) 2 + 0„03645A 4 (a+b) 4 - 0.00167 

A S (**b) 6 - $6 


(4.46) 



178 


Since a + p « 1.0. 

It Is to be noted that under steady state situation the 
generating solution becomes 

x 0 < t I,x) * A, a constant 

.Since | P | < 1*0 along the side AC of the stability triangle 

For this reason It is not necessary to collect the secular 
terms corresponding to (-P ) 71 term. 

Then y * A(a+b) with | y| < 3*0 

with the above constraint the eqn. ( 4 « 46 ) gives the following 
region 

0.98 < a - + b < 3*02 . 

Then the complete region In a-b plane In whiah Pj limit cycle 
oscillations are possible is obtained by Increasing D and is- ' 
given by 

a + b > 0.98 . (4.47) 

Similar analysis can be carried out to obtain regions 
corresponding to Pgt P 3 , P^ etc. limit cycle oscillations by 
suitably selecting a^ and ag in the basic system described in 
eqn. ( 4 . 36 ) so that the linear response is of periodic form 
with required, period. 

She various regions, showing P^ limit cycle oscillations 
are indicated In Fig. 4.6. 



179 


4.5 Simulation Results t 

Xhe theoretical findings given in the previous sect ion 
are verified by simulating the given system as a recurrence 
relation. Pig. 4*6 shows different regions in the a-b para- 
meter plane where limit oyole oscillation of different 
periods to exists. Figs* 4.7(a) - 4.7(f) show the time 
response for different set of [a,b] values taken from different 
regions shown in Fig* 4.6. Sots that Fig* 4.7(f) shows 
P 10 oscillation for a particular [a,b] values* This response 
changes considerably to another periodic response for small 
change in [a,b]* . 

The stability aspects of limit cycle osoillations are 
discussed in the following section* 

4.6 Stability of Limit Qy ole Osoillations i 

The simulated results obtained In the previous seat ion 
show that . the maximum Magnitude of the limit cycle oscillation 
is unity irrespective of its period of oscillation* This is 
true due to the saturation nonlinearity which damped the . 
output at + l.p whenever it exoeeds unity. The stability of 
limit cycle osoillations can be studied through variation . 
techniques. A small variation or perturbation is assumed around 
the known periodic limit cycle os dilation and a variational . 
equation is obtained* This equation is used to study the 
stability of limit cycles whose output sequence is known 
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beforehand. Ae mentioned in the introductory ohapter the 
stability of a limit cycle falls into any one of the following 
categories. 

(i) Stable limit cycle < A stable limit cycle is one for 
which the solution to variatlaval equation tends to 
zero as k increases to infinity. 

(ii) Semistable limit cycle : A limit cycle is known to be 
semistable if the variational equation exhibits dire- 
ctional stability properties. 

(lii) unstable limit oycle : The limit cycle oscillation is 
unstable if the variational equation gives an unbounded 
response. 

She variational equation for stability study of limit 
cycles In digital filters is obtained as follows s 

Let x*(k) be a known sequence representing a known 
limit cycle oscillation. Then introducing a small variation 
about x?(k) the filter variable x(k) oan'be expressed as 

x(k) ® xf*(k) + &x(k) 

I 

m 

where 6x(k) IS small perturbation about the limit cycle 
os cilia t ion. Substituting this In filter equation (4.1) 
we obtained. 

x*(ta*2) + &r<k+2) ■ f[ax*(k+l) + bx?f(k) + Ax(k+1) £ 6x(k)] . 
from which, for a positive valued perturbation 
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6x(k+2) - f[aa*(k*l) + bxi*<k) + adx(k*l) + bdx(k)] - a^(kf2) 

(4*48) 

and for a negative ‘valued perturbation 

Sx(kf2) - + bx*(k}- «&r(kfl) - b6x(k)] + x*(kf2). 

(4.49) 

The equations (4,48) and (4*43) can be ueed to study the 
. stability of given limit ay ole oscillation. As mentioned in 
the beginning, if 6x(k) * o. as k •** in both the equations 
the limit oyole is a stable limit cycle, if 6x(k) ♦ M, as 

k -*• , where H is finite, then the limit cycle is unstable, 
whereas if 6x(k) o as k** in one equation and 
fix(k) * K as k *■ in the other equation then the limit 
cycle is said to be semistable. 

For parameter values outside the stability triangle the 
linear response is unstable for all initial conditions. Shis 
kind of unstable response causes limit cycle oscillations in 
nonlinear case with saturation type of nonlinearity for almost 
all sets of [a,b] values,, except at oertain points which will 
be discussed later on in the analysis. The simulation result 
of the filter eqh. (4*1) shows that whenever there is a 
limit oyole os dilation, for the parameter values outside the . ... 
stability triangle, the limit oyele is always stable and 
' ^ this prope rty earn be oheoked by ocaslA sriag varicfue emamplee. 

For convenience and easy reference the F^. limit cycle 
oscillations given in Fig. 4*7 are considered for the stability 
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(i) Fig • 4.7(a) gives I., limit cycle with, following 

sequence 

x*(o) * 2.0 

I*(l) « 1.2 

J*(2) -**(5) - ...... -1.0 

considering the Initial disturbances 6x(o) * o.l v 6x(l) * 0.1 
the following values can be obtained for 6x(k) in the eqn. 
(4.48) 

6x(2) ■ f[3.6 + 4.04-0.3 + 0.2] - 1.0 * 0.0 
6x(3) - f (3. 04-2*4+0.0 + 0 . 2 ] - 1.0 = 0.0 

and 

6x(k) *0.0 Vk» 4. 

which shows 6x(k) ♦ o as k • 

and similarly from the eqn. (4.49) 

6x(2) « -f(3.6 + 4.0 - 0.3 - 0.2] + 1.0 - 0.0 
; dx(3) = -f [3. 04*2.4+0.0 - 0,2] + 1.6 *0.0 

fix(k) ■ 0 V" k > 4 

this also shows 6x(k) o as k . Therefore the 
Unit cycle is a stable limit cycle. 

(11) ?2 I*®** cycle is sketched in Fig. 4.7(b) and its 
sequence is 

x^(o) ** 2.0 
X*(l) ? 1.2 
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x*(2) - 1.0 

■ ■ j 

s*f3) - 0.4 
a^(4) - 1.0 
j^(5) - -1.0 

and x*(k) « 1.0 fop k even 

■. -1.0 for k. odd 

Then with 6x(o) » 0.1 and fix(l) = 0.1 and using eqn. (4.48) we 

obtain 

fix(2) » f[-2.4 + 4 - 0.2 + 0.2] - 1.0 *= 0.0 

&c(3) * f[-2 + 2.4 + 0.2] - 0.4 ■ 0.2 

6x(4) « f[-0.8 + 2.0 - 0.4] -1.0 « -0.2 

M5) - f[-2 + Q.8 + 0.4 + 0.4] + 1.0 a 0.6 

6x(6) - f[2.0 + 2.0 - 1.2 - 0.4] - 1.0=: 0.0 

■ _ ■ . . I 

fa(7) - f [-2.0 - 2.0 + 0.0 + 1.2] + 1.0 « 0.0 
6z(8) * 6x(9) * .... ■ 0.0 
that Is fix(k) * o as k +• 

and a similar result oan be obtained with eqn. (4. 49).. Xhls 
Shows the limit ay ole oscillation shown in Pig. 4.7(b) is stable: 

In a similar way the sequence 6x(k) for other periods of 
limit ay ales are obtained and plotted. In Figs. 4.8 and 4*9. 

Fig. 4.8 Is the time response of the variational equation 
(4.48) and Fig. 4.9 is that of .(4.49). .. The plots Show the 
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limit cycle oscillations are stable. 

The dependence of the limit cycle oscillations (for 
parameter values Inside the stability triangle ) on the Initial . 
condition has been briefly remarked by Olaasen et al [79] and 
Willson [91]. The limit cycles outside the stability triangle 
are also seen to be dependent . on the initial conditions. This 
shows that the tram 'limit cycles' in digital filters has a 
different meaning than that for nonlinear continuous time 
system. 

It Is observed that, due to the Instability behaviour of 
the linear response for the parameter values outside the 
stability triangle, the assumption that a limit cycle oscilla- 
tion is always possible Is not true with, saturation nonlinearity. 
For given parameter values [a, b] outside the stability triangle 
there exists a set of Initial conditions [x(o), x(l)J such that 
the nonlinear response dies down to origin monot onl oally . In 
otherwozds for particular filter coefficients outside the sta- 
bility triangle the linear response Is unstable whereas the 
nonl inear response with saturation nonlinearity is asymptoti- 
cally stable. This Interesting phenomenon is explained la the 
following sedtion. 

4.7 Monot onl o Decaying Besponsei Outside the Stability Triangle i 

. It is well known that the response to a linear second order 
difference equation Is asymptotically stable for the parameter 
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[a,b] values within the stability triangle, generally bounded 
when they are on the sides of the triangle and unstable for 
their location outside the triangle. Once the saturation 
nonlinearity is considered the response to a seoond order 

digital filter, for some [a,b] values located Inside the 

* 

stability triangle, be o ernes unstable and leads to self 
sustained oscillations [76, 79, 83]. This interesting inves- 
tigation tempts one to investigate the possible exist enoe 
of asymptotically stable response for the parameter values 
outside the stability triangle with saturation arithmetic. 

She details of the above investigations are given below s 


4.7.1 Analysis s. 


Tor a given a and.b, a set of initial values [x(o), 
x(l)] is obtained for which the response is asymptotiaally 
stable. Consider x(l) related to r(o] as 

x(l) ■ m x(o) 

where 


| m| < 1.0 


Then 


x(2) > f[ax(l) + bx(o)] 

«f[(anH-b) x(o)] 

How x( 2) can newer exoood unity In absolute valuo, that is 


|x(2)| » 1.0. for | (am+b) x(o) |> 1.0 

and x(2) ■ ■ (am + b) x(o) otherwise. . 
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z(3) * f[(a^n + ab + tm) z(o)] 

* [(a 2 + b) a + ab] x(o) 

because of the assumption of monotonio decay. 
z(4) ■ [ifi? + 2ab) m + a^ + b 2 ] x(o) 


x(N) = [(« 1 a 3r “^’ + a^a®"^b + b 2 + • • • •)m 

! 

+ + p 2 b 2 + Pja*"V + ...)]z(o 

- ■ i 

‘ " . ! 

Now, when. N is sufficiently large 

x(N) “ 0.0 

and the slope m can be evaluated from the right hand side of 
the equation for z(N) as i 

(4.50) 

where and p^ are the ooefflolente and ate given In table 
4.1 and 4.2 for different N. 

She elements In tables 4.1 and 4.2 are obtained as follows 
let 

-j (1 ■ 1,2,3, . ..) be the elements In the first . column. 
Then the Sable 4.1 can be filled -up with a i1 a l.oYland 


[P^"^ + Pga^b 2 + Pga N "^b 3 + ....] 

[a-j^a^"^* + «^a^“^b + gja*^b^”+ 777777] . 
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Sable 4.1 . 

Values of oc^ for different IT 
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Sable 4.2 

Values of for different ZT 
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and using the relationship t 

a JK-2(3-l),3 - E ^2(3-2), 3-J., IaiL * 2 »3» • •••* 
bd 

3 “ 2 » 3 | 4 » •• • • 
where 3 stands for column. 

Xa the same way elements of table 4.2 can also be obtained 
with a 11 = 0.0 and a^-, = 1,0 V i * 2 and 


a Bf2(J-i),3 



**■2(3-2), 3-1, 


M f 2,3| •••• 


3 * 2j3»4» •••«! 


. Before examining the detailed analysis, it is useful to 
outline the steps involved with the help of following flow abarts 


The flow dhart 4.1 explains the method of obtaining the 
initial condition locus whan the [a,b] values are specified, 
whereas the flolr chart 4*2 gives the Bteps to be followed to 
obtain a set of [a,b] values from a given [a 9 b] for monotonlo 
decaying kind of response. The set of initial conditions, for 
the different [a, b] values ore 3hst translation, of the initial 
condition lool obtained from flow chart 4.1. 

The value of n is known for a given [a,b] from the eqn. 
(4.50) 

Then, 

x(2) - (Sttfc) *(i) » (aart) i(o). 
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Flow Chart 4.1 Initial condition Looi for 

Honotonlo dooay respons e 






95 







196 


Nov the ma-rinmni amplitude that x(2) oan assume is 1.0 then 
3^(1) the value of x(l) oorree ponding to the maximum value of 
x(2) Is obtained as 

i - 

from vhloh 

V 1 * 88 amPE U»5D 

Similarly 

S <o) ” Swfc (4 * 52) 

again from 

x(2) - f[ax(l) + bx(o)] 

for ax(l) + bx(o) > 1.0 

x(2) - 1.0 • 

Under this condition x(l) » at^(l) » gg B g 

that is 

;g£g+tx(o) > 1.0 

* 

bx(o) > (4.53) 

am+b * 

likewise for 

ax(l) + *x(o) < -1-9 

. x(2) -1.0 

for this condition x(l) «° 3^(1) '« 
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then N*PO 

also whon jax(l) + bx(o)| 1 1.0. (4*54) 

x(2) * ax{l) + bx(o) • 

this gives the liter operation In the l i nea r range of the 
saturation nonl inearity , that Is 

*<2) - <***> *(D 

which Implies x(l) « nx(o) (4*5$) 

Case. 1 b positive 
from eqn. (4*53) 

*(°) - - V °) aai ■ z&b ( *; 56) 

. from eqn. (4.94) 

x <°> S aw* “ ~V o) “ 4 V 1 * * sSSs C4-.57)' 

The conditions given in eons* (4.55)# (4.56) and (4.57) are 

. : 

sketched in Fig. 4*10 for a positive and In Fig, 4*11 f or m 


*a (1) " wt (4.58) 

**(!> “kSb (4.59) 


negative. 

Qase 2 b negative 

from eqn. (4.53) 

*(°) - ssJe - *a(°) ana 
from eqn. (4.54) 




Xhe above conditions along with the condition given 
in eqn. (4.55) are plotted in Figs. 4.12 and 4.13 for m 
positive and negative respectively. 

Xhe Figs. 4*10 - 4*13 give the set of initial values 
[x(o) f x(l)] for a given [a,b] in the parameter plane for which 
the filter response decays to aero monot oni oally^ Xhe Initial 
condition loci in Figs. 4*10 - 4.13 give a mono tonic decaying 

response not only for a specified [a,b] value, but also for a 

■ ^ 

set of t*»b] values oonneoted with .slope m as given below s 
x(2) * f[ax(l) + bx(o) ] 
when the filter Is just nonlinear for 

x(2) « 1.0 » ax(l) + bx(o) 
or x(2) =-1.0 aaz(l) + bx(o) 
also we know for 

1( 2) -1.0, 1(1) - ^(1) - Sjjfe 

and for 

x(2) —1.0, *(1) — «^(1) "Safe 

then 

l.Q . (Stt) ^(1) 

from which 



r : :s$rjy 


5; v 
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Bqp* (4*60) gives a set of values [a,b] .for which there exist 

* 

a set of initial conditions [x(o), x(l)] as indicated in Figs. 
4.10 t 4*13 for which the nonlinear response of the filter 
decays down to origin monot oni cally . She region in the a-b 
parameter plane in which such type of. monotonic decaying res- 
ponse is possible is obtained as follows . * 

x(2) b ax(l) + bx(o) 
x(l) « mx(o) , 

Dividing one by the other 



then 



that is 


< 1 from the assumption 


< i.o 

® ' 

m nay be positive or negative, but | m| < 1.6. 
For m positive 


am+b < m 
b < (1-a) m 

•X < m < 1.0 



(4:61) 
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For m negative 

an + b < -a 
b < - (l+a)m 

-m < 1.0 . (4.62) 

Eqns. (4.61) and (4.62) give a region In a-b plane in which 
the solution deoays monotonically to the origin. This region 
is shown In Fig. 4.14. 

Note that the above results are not in conflict with 
earlier discussed limit cycles since only for a particular 
set of Initial •conditions determined by the Figs. 4.10 - 
4.13 will there arise a monotonic decaying sequence. In 

i 

general the family of limit cycles is much larger than that of 
the decaying sequence for a specified [a,bj. It is interesting 
to note that there exist regions outside the stability triangle 

where there cannot arise any monotonically decreasing response* 

■ 

■ " L 

She above discussed phenomena have been exactly verified 
by digital, computer simulation and few results. are given below. 

4.8 Examples t 

(i) For a ■ -a. 5 and b ■ 1.0 

Using eqn. (4.50), the value of m is obtained as 0.5. 

■* " . 

For this {a y b] t the set of Initial conditions that are 

M ■ " " _ B _ , 

" I . B ■ 

responsible for monptonio decaying response is given in 
Fig. 4.13. Whereas.. the Figs* 4.16(a) - 4.16(d) show the time 
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response of the system for different Initial values [x(o), 
x(l)] taken from Pig. 4-15* that is 

Pig. 4.16(a) , for x(l) = -2.0 and z(o) - -4.0 

Pig. 4.16(h) , for z(l) = 2.0 and x(o) - 4.0 : 

Pig. 4.16(e), for z(l) * -1.5 and z(o) - -3.0 

Pig. 4.16(d); for z(l) - 1.0 and z(o) * 2.0 

(11) Por a * -2..0 and h '■ -0.75 

the value of m Is obtained as -0.5 

for this [a,b] value the Initial condition locus is 

■■ 

sketched In Pig. 4.17. 

The tine response for various ohosen Initial conditions 

i ■' ■ 

from the Pig. 4.17 is shown In Fig.. 4.18 v that is 
Pig. 4.18(a) , x(o) > -4.0 and x(l) = 2.0 

* 1 

Pig. 4.18(b), z(o) ^ 4.0 and x(l) =-2.0 

Pig. 4.18(c), x(o) = 2.0 and z(l) =-1.0 

Pig. 4.18(d), x(o) = -2.0 and x(l) = 1.0 

Thus the theoretical findings are exactly verified with 
simulated results of the given system. 

4.9 Conclusion i 

B ■ I , - 

The existence of limit cyole oscillation In a second 
order digital filter, vith overflow nonlinearity of saturation 
kind, Is demonstrated with the parameter values outside the 
.stability triangle In the a-b parameter space. Three methods 
have been proposed for. locating regions in the parameter plane 
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FIG- 6-15 INITIAL .CONDITION LOCUS FOR cu45,b=1-0 AND m = 0-S 
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where there exists limit 07 ole oscillation of a specified 
period. Xhe first of these is a variation of one proposed by 
Ebert et al [76] and the other two are based on a polynomial 
approximation to the saturation nonlinearity. In these two 
methods a lower degree polynomial with small range for the 
dependent variable may be considered and then by induction one 
can obtain the complete region where limit cycle oscillation of 
specified period is possible. Xhe third method exploits the 
use of multiple scale perturbatlonal approach and is. believed 
to be quite novel. 

Regions in the parameter plane where the output sequence 
decays monot oni oally to zero have also been ascertained. 

These are functions of the initial conditions and an explicit 
relation linking the parameter [a,b] and the initial data 
has been determined for this, purpose. These relations both 
in the parameter space and in the initial condition spaoe 
have been appropriately, mapped. 



CHAPTER 5 


JUMP PHSNOMBMQM US’ DIGITAL FILTERS 
5*1 Introduction i 

She proceeding chapter has considered In detail the 
analysis of seoond order digital filters with saturation 
overflow nonlinearity under force free situation. The 
overflow nonlinear phenomenon Introduced in the normal 
operations of the filter gives rise to limit cycle oscilla- 
tions of different periods for parameter values outside the 
stability triangle. She stability analysis of such limit 
cycles was carried out using a. variational technique. Fur- 
ther investigation had also revealed the existence of a 
region in the a-b parameter spaoe outside the stability tri- 
angle where the output sequence decays down to zero mono- 
tonloally with time. In this chapter a seoond order digital 
filter with an external input is considered. If an input 
signal is applied to the filter then the situation is 
different and more complicated than for the foroe free 
situation. For periodic input signals the quantization errors 
will be periodio too, and are sometimes referred to as limit 
cycles [93] and the nature of such limit cycles will be 
highly dependent oh the specific input signal and the quan- 
tization effects are therefore not considered in this ana- 
lysis; As mentioned earlier, the overflow though occur lng 
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rarely in a proparly designed digital filters Introduces a 
large error In the normal operation of the filter and are the 
therefore Investigated farther in this chapter. The study, 
of overflow oscillations m the forced response Of a second 
order digital filter was initially' due to Willson [94], who 
shoved .that unlike the zero input case, the response, for sobs 
choice of filter coefficients [a,b], differs from the ideal 
linear response and the filter turns out to he asymptotically 
unstable. She region in the a-b plane where this instability 
arises is shown in Fig. 2 of [94]. This result has recently 
been significantly extended by Clausen and Eristiansson [112], 
The following section is concerned with the jump phenomena in 

a f oroed second order digital filter with overflow saturation 

% 

nonlinearity. This jump phenomenon will be shown to be closely 
related , to the overflov nonlinear oscillation described in 


[ 94 , 112 ]. ’ 

9.2 Jump Phenomenon . in Digital Filters i 


A particular kind of non! inear effect known ae the jump 
phenomenon has been examined recently by Krlstianseon [95] 
and Claasen and Eristiansson [96]. A somewhat related study 
is due to Willson [91] who has obtained different regions 
inside the stability triangle where foroed overflow oscilla- 
tions of different periods can be sustained. There appears 
to be no other speolflo study of this jump phenomenon in 


.IgUalje 


liters in the available literature. 
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Jump phenomena are characterised by changes- in the 
steady state amplitude of the filter response that may occur 
due to small change In the amplitude or frequency of the peri- 
odic input signal* The present . study is. concerned vith further 
investigations of the jump phenomena In a second order digital . 
filter vith a saturation nonlinearity arising due to overflow 
in the adder* She study provides information about the range 
of the input amplitude for which there exists a jump* For 
given values of the filter coefficients [ a t b] and input 
frequency a, conditions are derived which guarantee a jump or 
sustained nonlinear oscillations* A graphical construction of 
the derived conditions provides the location of the regions 
inside the stability triangle where such jump behaviour and 
nonlinear sustained oscillations are possible* 


5*3 Sjyeitem Model t 


The second order section of the digital filter under 

study is given In Fig. 5*1 and is described lay 

\ 

x(k+2) * f (ax(k+l) + bx(k) + s(k)) (5*1) 

where f ( •) is the nonlinear function defining the 
saturation arithmetic by 


f(a) « 


e for | a | < 1*0 
sgn(c) for |.o t > 1*0 



and the graphical representation is shewn In 7ig. 5«2. Z(k) 
Is the external Input to the filter, generally periodic and 
is given by 

Z(k) ■ 7 oob ok 

where 

7 is the input amplitude and a is the input frequency. 

Of Interest is the examination of conditions under which 
disturbances (equivalent to small changes) in amplitude of 
input at specified times lead to a jump from one steady state 
(assumed periodic) solution to another steady state solution 
hut with a different amplitude • Such an analysis is carried 
out in the following section. 

5.4 Analysis : 

Xhe linear version of the nonlinear system described in 
eqn. (5.1) is 

x(k+2) = ax(k+l) + bx(k) + Z(k) ($.2) 

whose solution is 

x(k) » x^(k) + Xg(k) (5.5). 

where 

x^(kj is the transient solution obtained assuming sero 
input 

x g (k) is the contribution due to the forcing function. 




FIG- 5-1 FORCED SECOND ORDER DIGITAL FILTER 
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As .pointed out earlier, the Interest 
is to study the iwaar. effects in a second order digital 
filter with saturation overflow nonlinearity with external 
forcing function for the coefficient values Inside the 
stability triangle. This implies the transient response 
Xj.(k) * o aa k «» irrespective of Initial conditions and 
XgOO will be the steady state solution of eqn. (5.2). 

Xhe transient term Xj.(k) In eqn. (5.3) is obtained ae 
follows t 

let z(k) = 0, then the system equation (5-2) takes the form 

x(k+2) - ax(ktl) - bx(k) = 0. 

3?or this equation, 

X^-ax-baO 

is. the characteristic equation, from which the roots are 
obtained as 

x i,2 "! ± ^<i> 2 + b - 

She roots are oscillatory (complex conjuoqte) type if and 

tL 2 

only if b is negative ...and | b | > (§} , that Is 

X i#2 ■ | ± 3 l/ -b - (|) 2 
from which 

'HI ■ r » f-fc 
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and cos p « , tt»* 1» a - 2r oob 0. 

where r and p are the amplitude and phase of the root* 

Thus the transient solution is given by 

x t (k) m x r k oos pk + bJ* sin pk , (5*4) 

Here A and B are the constants to be evaluated from the initial 
data x(o) and x(l). 

The steady state term Xg(k) in (5.3) Is evaluated as 
follows i 

Let . s(k) * B cos ok . 

Assume the steady state solution of (5.2) as 

x 8 (k) = C oos ok + D sin ok (5.5) 

where 0 and B are the constants that are related to the Input 
variables namely* the amplitude B and the frequency a and 
their determination proceeds as follows * 

Bram eqn. (5.5) 

Xgtktl) * (0 oos a + D sin a) oos ok - 

(0 sin a - B oos a) sin ok 

x B (k*2) ■ (C oos 2s + B sin 2a) oos ok - 

(C sin 2a - B oos. 2a) sin ak , 
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Substituting the above expressions in eqn. (5.2) and 
ooUeotlng the like terms, the following equations are 
obtained* 

Cfcoe 2a - a oos a - b] + D[stn 2a - a sin a] - 7 * 0.0 
0[sin 2a - a sin a] * 2)[oos 2a - a cos a - b] *» 0.0 
from which 

0 = F(oos 2a - a oos g - h) 

1 + a Z +b Z - 2a(l-b) oos a - 2b oos 2a 

p aiflug)., 

1 + a z +'b z - 2a(l-b) oos a - 2b oos 2a * 

m ■ 

How the complete solution is 

* 

x(k) * Ar* oos pk + Br^ sin pk + 0 oos ok + D sin ak« 

The transient part of the above solution can be completely 
eliminated by proper choice of the initial conditions x(o) and 
x(l), since 

x(o) « A + 0 

x(l) = Ar oos p + Br sin p + 0 cos a + D sin a 
and by ahoosing 

x(o) ■ 0 » P{oofl 2a - a oos a - b)/fi^ (5*6) 

x(l) * 0 oos a + B sin a = F[ (1-b) oos. a - aj/B^ (5*7) 
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■where 

B^ « 1 + + b® - 2a(l-b) cos a - 2b cos 2a. (5*8) 

z(k) asm now be rewritten as 


x(k) a o eoe ode + D sin ode 


£ oca (ode + 0) 


where 


E 


5 


« - tan -1 fS** Iff. - »■ Sfrff-g.. i 

L 008 2a - a 008 a-t> J 


t5.9) 


y ■ * 

The following renark oan now be made. In order to 
cause a jump from a linear oscillation (satisfying f(e) = a) 
to a nonlinear os oil la t ion (where the saturation ohazact eristic 
is effective) It is necessary that 


a < i.o (5.10) 

and E being real constrains B^ > 0. 

Further, If half a period of the input sinusoid is 
assumed to be an Integer times the sample interval [951 then 

^ si (Integer XT). 

Consider now the search for a particular type of nonlinear 
oscillation which has same frequency as Input frequency and 
will remain damped at the saturation level + 1.0 at four 
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sampling instants per ay ole (first two successive instants 
during positive half cycle and other two successive Instants 
during negative half cycle). A mathematical statement of the 
prescribed nonlinear oscillation is 

x(k) *1.0 k * 0,1 (5.U) 

x(k) * ax(k-l) + bx(k-2) + 3? cos <x(k-2) k = 2,3, ...ft-1 

( 5 . 12 ) 

x(k) * -1.0 k * IT, ft+1 (5.13) 

She above description essentially implies that the 
system is oscillating in the linear mode (l.e. R < 1.0) 

before disturbance application. She disturbance is injected 

3 

at say k * 0,1 oausing the response x(k) to have unit value 
for k * 0,1. 

It will now be shown how a particular choice of ft 
enables determination in the a-b parameter plane, regions where 
transition from a linear solution with R < 1.0 to the above 
mentioned nonlinear oscillation take place. It is apparent 
that combination . of conditions given in eqns. (5.11) to (5.13) 
together, with the understanding that R < 1.0 from eqn. 

(5.10), should provide the conditions for jump to ooour* 

She method of deriving the required conditions is given 


below for various ft 
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Period 6 pagination (H * 3) s 

For til * 3* the foil owing conditions are readily 
% - 1 + a 2 + b 2 - a(l-b) + b > 0.0 

x(k) - 1.0 k * 0,1 

x(k) ** ax(k-l) + bx(k-2) + F cos y (k-2) k = 
x(k) » -1.0 k * 3,4 

. Thus 

x(o) *1.0 
x(l) * 1.0 
x(2) * a + b + F 
x(3) * -1.0 
x(4) = -J..0 

How from eqn. (5.1) 

x(3) * f[ax(2) + bx(l) + F oos */3] 

Thus using the earlier derived conditions 

-1.0 = f[a(a+b+F) + b + F/2] 

which Implies 

a(a+b+P) + b + F/2 < -1.0 

from whioh. 

(a+0.5) F < - [a(a+b) + b + 1.0] # 


obtained 


2 

( 5 . 34 ) 


( 5 . 15 ) 


( 5 . 16 ) 
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Further 

z(4) - f[az(3) + bz(2) + F oos 2*/3] 

m 

that is 

' -1.0 * f[*4 + 1) (a+b+F ) — F/2] 
from which 

-a + b (a+b+F) - F/2 < -1,0 

(b - 0.5)F < - [b(a+b) - a + 1.0] « (5.17) 

i 

Thus we can state the necessary and sufficient conditions for 
a jump from a linear steady state oscillation with magnitude 
lees than unity to a nonlinear os dilation (described above) as 


0 - 1 t 

F 2 < ^ 


0-2 f 

(a+0.5)F< - {a(a+b) + b + 1.0] 


0 - 3 I 

(b-0.5)F < - [b(a+b) - a + 1.0] 

(5.18) 


It is also obvious that satisfaction of conditions 0-2, 

0-3 alone provides only the condition for sustained os dilations 
without necessarily the presence of a jump. Conditions 0-1 
provides the additional desired jump condition. 

Period 8 oscillation (K - 4 ) s 

Vlth N *» 4, the assumed relations in (5.11) - (5.13) 
assume the following form 
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x(k) -1.0 k « 0,1 

x(k) - ax(k-l) + bx(k-2) + 7 oos J (k-2) k - 2.3 
x(k) - -1.0 k — 4*5 (5.19) 

Then 

x(o) = 1.0 
x(l) - 1.0 
x(2) - a + b + 7 

x(3) » a(a+b+7) + b + 7/1T2 = a(a+b) + b + (a+i-)7 

. V2 

x(4) - -1.0 
x(5) - -1.0 

from the system equation (5.1) 

X(4) - f[ax(3) + bx(2)] 

from the above assumed relation x(4) - -1.0, we obtain 
-1.0 - £[ax(3) + bx(2)] 
that is 

ax(3) + bx(2) < - 1.0 

which implies 

a 2 (a+b) + ab + a (a + y^-)7 + b(a-fb) + b7 < -1.0 

7 (a 2 + + b) < — (a 2 (a+b) + 2ab + b 2 + 1.0]. (5.20) 

Further 

*(5) - - 1.0 - f[ax(4) + bx(3) - ?/fa] 
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that is 

ax(i) + bx(3) -F/Y2 < -1.0 

-a + ab(a+b) + b 2 + b(a + - V/f2 < -1.0 

front which 

F(ab + b/Y2 - < - [ab(a+b) + b 2 - a + 1.0] (5.21) 

Then the necessary and sufficient condition for a junp 
to take place is 

0-* I» 2< \ 

0 - 2 t J(a 2 + p + b) < - [a 2 (a+b) + 2ab + D 2 + 1.0] 

0 • 3 ) ^ . * [ai(a+b) + 1> 2 *. a + X«0 j 

(5.22) 

where 

^ - 1 + a 2 + b 2 - V2 a(l-b). 

As in the ease of period 6 oscillations, conditions 
0-1, C-2 and 0-3 are the necessary and sufflolent conditions 
for jump to take place and 0-2 and 0-3 are the necessary and 
sufflolent conditions for the nonlinear oscillations to 
sustain. 

Period 10 oscillation (S«5) t 

Proceeding In the sane way as above the following 
relations are obtained ' 
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x(k) • 1.0 k - 0,1 

x(k) » ax(k-l) + bx(k-2) + 7 oos g(k-2) k - 2,3,4 
. x(k) - -i,0 k-5,6 <5.23) 

from which 

X(o) m 1*0 
X(l) « 1*0 
x(2) *= a + b + F 
x(3) ■ a(a*b) + b + P{a + 0.809) 
x(4) * (a 2 +b)(a+b) + ab + P(a 2 + 0.809a + b + 0.309) 
x(5) - f[(a 3 +2ab)(a+b) + a 2 b + b 2 + F(a 3 +0.809a 2 + 2ab 
+ 0.309a + 0.809b - 0.309)] 

x(6) - f[(a 2 b»>b 2 )(ari>b) *i«b 2 - a ♦ Ffa^OrfSOOab + b 2 

■ ■ 

+ 0.309ab r 0.809)] # 

Then imposing the oondltiona aa x(5) and x(6) given, la (5.23). 
ve have 

[ (a 3 +2ab)(a+b) +■ aPb + b 2 + F(a 3 + 0.809a 2 + 2ab + 0.309a 

+ 0.809b - 0.309)] < -1.0. 

[(a 2 b+b 2 )(a+b) + ab 2 - a + ffah + 0.809ab * b 2 + 0w?09b 

- 0.809)] < -J-.0. 

She necessary and sufficient oondltloas for jump phoneme** 
non nonlinear sustained oscillations can be obtained as . 
follows from the above equations. 
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0-1 s F 2 < 

0-2 i (a 3 + 0.809a 2 + 2ab + 0.509a + 0.809b - 0.509)? 

< - [1.0 + fe 3 + 2ab)(a+b) + a^ + b 2 ] 

0-5 : (aS + 0.809ab + b 2 + 0.509b - 0.809)i 

< - [1.0 + (a^b + b 2 )(a+b) + ab 2 - a] 

where 

\ « 1 + aib 2 - 2a(l-b) eos */5 - 2b cos 2a/5. (5*24) 

Period 12 oscillations (H*=6) : 

She conditions stated in eqps. (5.11) - (5.15) take the 
following form for 9*6. 

ar(k) -1.0 k - 1,2 

x(k) - ax(k-l) + bx(k-2) + P cos g(k-2) k - 2,5,4,5 
x(k) - -1.0 k - 6,7 (5.25) 

She above procedure oan now be followed to obtain the 
conditions 0-1, 0-2 and 0-5. 

0-1 i P 2 <■ 

0-2 r (a 4 + a 2 + Sa^ + §- + Y5ab + b 2 + \ - *)F 

« - [1.0 + (a+b) (a 4 +5a 2 bfb 2 ) + ab(a 2 +2b)} 

0-5 . i (a^b + ^ a.h» + 2ab 2 + S| + *Jb 2 - 

.< - [1.0*(a+b) (a 3 b+2ab 2 ) + a^^b 3 -*] 

(5.a) 



where 

- 1 + a 2 + b 2 - i(3a(l-b) - b. 

Here again all the above three conditions are necessary 
and sufficient to cause a jump and the last two conditions 
are necessary in order to cause nonlinear sustained os dilations . 

i 

!Che conditions stated In eqns. (5.11) - (5-15) are the most ; 
general for all N. For ff ■ 3,4,5 and 6, the eqps» (5.18), (5.22) 
(5.24) and (5.26) respectively provide the necessary conditions . 
From these conditions the following useful results oan be 
deduced* 

(1) For given input variables (F and a), the values of filter 
coefficients [a # b] satisfying all the three inequalities 
0-1 to 0-3, give a region in the a-b parameter plane in which 
a jump behaviour is observed. 

(li) For given input frequency and amplitude, the satisfaction 
of, conditions 0-2 and 0-3 alone gives a region where there 
exists a nonlinear sustained oscillation. 

(ill). For a given input frequency and given [a,b ] values, 
the values of F, if any, satisfying all the three inequalities 
0-1 to C-3 provide the range of F for which there exists always 
a jump, in the solution. 

She conditions 0-1 to 0-3 (assuming equality in their 
respective definitions) are plotted in the a-b plane far 
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various 'values of the Input amplitude Pig. 5-5 gives the 
oomplete plot of all the three conditions stated In (5*18) 
for V * 3# whereas Figs. 5-4» 5-5 and 5*6 provide the plot 
of all the conditions given In eqps. (5*22) , (5.24) and (5*26), 
all for a particular F. Zhe oomplete oomputed values of eqns» 
(5.22), (5.24) and (5*26) are tabulated In Sables 5*1*. 5*2 and 
5*3 respectively. 

For Illustration, for ET = 3 and for F * -0.5* the region 

ABC3DB in Fig. 5*3 provides a region where .the nonlinear phe- 

• . 

nomena exist. Zhe region GDtIO gives the region where the jump 
phenomenon is possible and the region ABGBA gives the region 
for nonlinear sustained oscillation. For other values of I, 
the above regions are indicated in the respective figures. 

Zhe line PQ in Fig. 5*3 and the curve PQ in Fig. 5*4 define 
the locus of the tip of the boundary of the region where the 
nonlinear sustained oscillation is possible. 

Once the values of input amplitude and the input frequency 

are known, the region in which, the jump phenomenon takes place 

■ ■ ■ 

can be determined as mentioned in. the previous section. Sere 
the jump behaviour is observed by applying the perturbation in 

" j 

the system variable euoh that x(j) ■ 1.0 for j * 2 Hn 9 
2Nn + i for ia particular integer n > 1. It is to be noted 
that the conditions 0-1 to, 0-3 are derived baaed on the assump- 
tion that x(o) * x(l) ■ 1.0 when, the nonlinearity is effective 
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Sable 5.1o CGHPUTHD 7AIUBS PQR 0QSDII10V 0-3* 1-4 


a 



7AHJBS of b 

P—0.4 7—0.3 

7—0.2 

7—0.1 ■ 

1.00 






-0.1200 

* 

! 





■ 

-0.2947 

1.05 






-0.4285 

1.10 





-0.3473 

-0.5171 

1.15 




-0.2309 

-0,4512 

-0.59U 

1.20 




-0.3800 

-0.5313 

-0.6573 

1.25 



-0.2975 

-0.4694 

-0.6008 

-0.7184 

1.30 



-0.4042 

-0.5433 

-0.6641 

-0.7760 

1.35 


„ 

-0.4844 

-0.6093 

-0.7232 

-0.8310 

1.40 


-0.4235 

*0.5537 

-0.6701 

-0.7793 

-0.8841 

1.45 


*0.4972 

~0«6167 

-0.7275 

-0.8333 

-0.9356 

1.50 

-0.5029 

-0.5627 

-0.6755 

-0.7884 

-0.8654 

-0.9859 

1.55 

—0.5668 

-0.6233 

-0.7314 

-0.8353 

-0.9323 


1.60 

-0.6263 

-0.6803 

-0.7851 

-0.8867 

-0.9860 


1.65 

—0.6826 

-0.7349 

-0.8371 

-0.9369 




-0.7365 

-0.7876 

-0.8878 

-0.9862 



1.75 

-0.7887 

-0.8388 

-0.9375 




1.80 

-0.8396 

-0.8889 

-0.9863 




1.85 

-0.8894 

-0.9380 





1.90 

-0.9383 

-0.9864 





1.95 

-0.9864 






2.00 
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Sable 5.2a 8ZMDB.S1D VALUES TOR OOHLITIOH 0-1, H « 5 


vEsx ns n j F T 


a ' 

F— 0.5 

T— 0.4 

F-0.3 

E— 0.2 

E-O.l 

1.00 

-0. 8**36 
-0.1564 

-0.6674 

-0.3326 




1.05 

-0.9126 

-0.1683 

-0.7607 

-0.5202 




1.10 

-0.9775 

-0.1843 

-0.8403 

-0.3215 

-0.5809 



1.15 

-0.2058 

-0.9117 

-0.3310 

-0.7410 

-0.5017 



1.20 

-0.2263 

-0.9774 • 
-0.3462 

-0.8539 

-0.4897 




-0.2515 

-0.3658 

-0.9101 

-0.4944 

-0.7023 

Vi 

1.50 


-0.4155 

-0.5279 

-0.9064 

-0.6600 

■ 

1.40 


-0.4447 

-0.5524 

-0,9772 

-0.6701 

-0.8236 

1.45 


-0.4764 

-0.5808 

-0.6902 

-0.8484 

-0.8797 

1.50 


-0.5106 

-0.6126 

-0.7169 

-0.8325 

-0.9765 

1.55 



-0.6476 

-0.7490 

-0.8533 

1.60 



—0.6856 

-0.7857 

-0.8860 

1.65 



-0.7265 

-0.8268 

-0.9276 

1.70 



-0.7702 

-0,8722 

-0.97867 

1.75 



*0.8170 

-0.9224 


1.80 



-0.8669 

-0,9782 

• 

1.85 



-0.9204 



1.90 



-0.9781 



1.95 




. 


2.00 

■ 


' 





Sable 5.2b SIMULATED VAKJBS OF CONDITION 0-2, S ■ 5 
| VALUES OF b 

■MM iMBMMMMMMaMMaMaMMMMMMWMMHaMHMMMMMlMMMMMMWaMMaMMaeHMaM 

gp-0.7 E«-0.6 P—0.5 F«-0.4 E^O.3 P«0.2 3M).l 


1.00 

1*05 

1.10 

1.15 

1.20 

1.25 

1.30 

1.35 

1.40 

1.45 

1.50 

1.55 

1.60 

1.65 

1.70 

1.75 

1.60 

1.85 

1.90 

1.95 

2.00 


-0.4291 

-0.6891 

-0.3785 

-0.8570 

-0.3625 

-0.9932 

-0.3604 


-0.5109 

- 0.8192 

-0.4717 

-0.9788 

-0.4601 

-0.4609 

-0.4693 


- 0.5892 

- 0.9560 

- 0.5616 

-0.5553 

-0.5595 

-0.5704 

- 0.5862 

- 0.6056 


-0.7336 

- 0.9062 

—0.6662 

-0.6494 

-0.6484 

-0.6564 

-0.6703 

-0.6885 

- 0.7100 

-0.7343 

- 0.7608 


-0.7774 

-0.7431 

-0.7358 

-0.7404 

-0.7522 

-0.7692 

-0.7899 

-0.8137 

- 0.8400 

-0.8684 

-0.8985 

-0.9302 

-0.9633 


-0.9124 

-0.8363 

- 0.8206 

-0.8216 

-0.8315 

-0.8473 

-0.8674- 

-0.8909“ 

-0.9170 

-0.9454 

-0.9756 


-0.9285 

-0>9020 

-0*8992 

-0.9073 

- 6.9222 

-0.9420 

-0*9653 

-0.9915 
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> Sable 5.2o SBRJIATED VAKJBS POE UK? CCHDITIOH 0-3, H -■ 5 
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Sable 5.3a SIMDXAISD 7ALDBS POE OOHDIIIQN 0-1, » * 6 





Sable 5.5b SXMDIATED VAIDES 


YAHJES 



1.15 

1.20 

1.25 

1.50 

1.55 

1.40 

1.45 

1.50 

1.55 

1.60 

1.65 

1.70 

1.75 

1.80 

1.85 

1.90 

1.95 

2.00 


-0.5000 

-0.5451 

-0.5905 

-0.4568 

-0.4859 

-0.5519 

- 0.5808 

- 0.6506 

-0.5412 

-0.5502 


-0.5588 

-0.5856 

-0.4290 

-0.4755 

-0.5224 

-0.5705 

-0.6195 

-0.6695 

-0.7205 

— 0.6561 

-0.6585 

-0.6444 

-0.6590 

-0.6791 

-0.7050 


6 6c 
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OCHDITIOIT 

b 




-0.5770 
r .4217 
.4672 
.5155 
-0.5607 
-0J6089 
- 0.6580 
-0.7081 
-0.7592 
-0.8115 
-0.7641 
-0.7498 
-0.7564 
-0.7721 
-0.7956 
- 0.8190 
-0.8474 
- 0.8782 
-0.9110 
-0 


0-2, U = 6 


F»-0.1 

-0.4147 

-0.4595 

-0.5050 

-0.5514 

-0.5987 

-0.6470 

-0.6965 

-0.7465 

-0.7978 

- 0.8500 

-0.9055 

- 0.8870 

- 0.8650 

-0.8715 

-0.8880 

- 0.9106 

-0.9575 

-0.9671 

-0.9992 
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Sable 5.3o SIMULATED VALUES FOR COHDITION C-3, E * 6 


A 


VALUES 

OF b 


P »- 0.3 


p ::.- 0 a 

1.00 


- 0.0582 

- 3.1113 

1.05 


- 0.0859 

- 0.2290 

1.10 


- 0.1121 

- 0.3392 

1.15 


- 0.2501 

- 0.4149 

1.20 

■ 

- 0.3589 

- 0.4764 

1.25 

- 0.25914 

- 0.4056 

- 0.5296 

1.30 

- 0.3352 

- 0.4619 

- 0.5771 

1.35 

- 0.3959 

- 0.5120 

- 0.6208 

1.40 

- 0.4488 

- 0.5580 

— 0.6620 

1.45 

- 0.4969 

- 0.6011 

- 0.7013 

1.50 

- 0.5419 

- 0.6423 

- 0.7395 

1.55 

- 0.5847 

- 0.6823 

- 0.7770 

1.60 

- 0.6261 

- 0.7214 

- 0.8142 

1.65 

- 0.6666 

■ — 0.7600 

- 0.8513 

1.70 

- 0.7064 

- 0.7984 

- 0.8885 

1.75 


- 0.8368 

- 0.9259 

1.80 


- 0.8752 

- 0.9636 

1.85 


- 0.9139 


1.90 


- 0 . 952 ? 


1.95 


- 0 . 99a 


2.00 


' 
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112]. Depending on the time Instants at whiah x( j) ■ 1,0 
.nd x(J+l) *» 1,0 are applied for j *» 2Kn , sometimes transients 
nay be observed when the response jumps from one steady state 
solution to another. In the present study the magnitude of 
the disturbance te be applied may be large and it Is. of. 
Interest to ascertain how small a disturbance can cause the 
above prescribed nonlinear oscillation. 

Minimum Disturbance required to cause nonlinear oscillations i 

Instead of disturbing the system response at time in- . 
starts j, j+1 suoh that x(j) = 1.0 and x(j+l) * 1,0, one can 
try to obtain the smallest values of the disturbances for 
which, the same response as above can be obtained under steady 
state. The procedure to obtain these conditions ‘are as 
follows I 

Consider the given system and let the initial values . 
be x(o) - x Q and x(l) = x^, which are to be determined such 
that the output sequenoe assumes the form given In eqns^ 

(5*11) to (5.13). Then, 

x{2) ■ axj_ + bx Q + P 

m 

x(3) * ax(2) + bx^ + P cos a 
x(4) ® ax(3) *1* bx(2) .+ P oos 2a 


(5.27) 
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Bqpu (5*27) is a sat of egressions for which j x(j) | 
is oonstrainsd to be less than or equal to unity. (Chat Is 
although each of the equations In (5. 27) describee a linear 
system, a properly combined version ;vay lead to a nonlinear 
oscillation. Then knowing the reglo:. In the a-b plane for .a 
specified Input amplitude and frequency, the required minimum 
values for x Q and x^ can be aalculate-d using the relations 
given in (5*27). For example. 

(1) Let H « 3# the period 6 oscillation 


She expressions In the eqn. (5.27) take the following 

form 

x(2) « ax^ + bx Q + F (5.2B) 

x(3) * -JL.Q .» ax(2) + bx^ + F cos w/3 « ax(2) + ba^ + F / 

(5.29) 


x(4) • -1.0 - ax(3) + bx(2) - F/2 = 

from the last expression 

x(2) - (a + F/2 - 1.0)/b 
and from equations (5.29) and. (5.28) 

‘ Xj/ ■ - (1.0 + ax{ 2) + F/2)/b 

x 0 « (x(2) - ax^ - F)/b 


+ bx{2) - F/2 

(5.30) 


(5.31] 

(5.32] 


She eqns. (5.31) and (5.32) provide the minimum values tn 
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and for XT ■ 6 

x(5) - (a + P - 1.0)/b 
x(4) - - (1.0 + ai(5) - P/2)/b 
x(3) - (z(5) - ax(4)/b 
x(2) - (x(4) - ax(3) - P/2)/b 

ij. - (x(3) - ax(2) - P)/b (5.38) 

x„ - (x(2) - axi - P)/b. (5.39) 

Then with these x Q and values and for the known filter 
coefficients [a,b] In the region In which the lump or non- 
linear sustained oscillations are possible, the time response 
of the system equation (5*1) is computed and the results are 
discussed in the following section. 

5.5 Simulation Results j 

The theoretical investigations reported in the above 
section are verified by actually simulating the system 
equation Sas a second order recurrence relation. The follow- 
ing are the results corresponding to the observations made in 
the previous section. 

5.5.1 Jump up phenomena s 

The filter parameters [a,b] are chosen from the region 
where the jumps are possible for specified input amplitude, 
and the minimum values of and required to cause the 
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nonlinear os dilations are obtained for various H as indioat ed 
in section 5.4 and using these values such that x( j) =* x Q and 
x( J+l) - x^ the time response of the given system is simulated 
for various discrete time index k. For all these oases the 
linear response, whose steady state amplitude is less than unity 
is evaluated using the Initial oonditlons x(o) and x(l) given 
in eqms. (5*6) and (5*7)< 

Fig. 5.3 shows a region GSDSO obtained in the a-b parameter 
plane, the interior of vhioh satisfies inequalities 0-1 to 
0-3 for 5=3. Ibis above region is obtained for F = -0.3. 

Then consider a » 1*6 and b = -0.9 inside ODBC, with these 
values and from eqps. (3.31) and (9.32), we have 

x Q « 0.152 

x^ « 0.182 (5.40) 

. and from eqns. (5.6) and (5.7) the following initial values are 
obtained 

x(o) • 0.43 

x(l) « 0.819 (5.41) 

Fig. 5.7a gives the time response with Initial conditions 
given In (5.41) and Fig. 5.7b shows the time plot with x(o) = x Q 
. and x(l) = x 1 , where and are given ln eqn. (5.40). Figs. 
5.7o to 5.7h show a jump up phenomena between the solutions 
when the x Q and values are applied suoh that x( j) * x G and 



» 
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e 


x(K) 


(e) Jump between solutions x (501-0-152, x (51 ) = 0-182 


x IK) 



(f ) Jump between solutions x (51)= 0-1 52, x (52) =0-182 


x(K) 



(g) Jump between solutions x (52)= 0-152, X (531= 0-182 


x(K) 


(h) Jump, between solutions x ( 53 )= 0-152 , x (5 2,1=0-182 
FIG. 5-7 SOLUTION TO EQN. (51) WITH q= 1-6 , b^O-98,F-.-0-5 
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x( j+1) - where j is the integer valued time instant. It 

la observed by simulation that for j ■ 2 UET (n = 1,2, 
a jump up phenomenon is observed and the steady state non* 
linear osoillations has the description given earlier after 
a few ay ole of transient behaviour. The Figs. 5.7o to 
5.7h show the jump behaviour between the steady state solu- 
tions when x Q and x^ are applied at various time instants such 
that x(j) » x Q and x(j+l) * for j = 2nXT + m, that is for . 
n = 6 and m » 0,1,2, .... 5 respectively. It is interesting 
to observe the transient behaviour that is developed before 
the steady state nonlinear osolllatlons sets in. 

For X « 4f similar results are plotted in Figs. 5.8a to 
5.0h with a = 1.85, b = -0.95, x(o) * 0.606, 'x(l) ■ 0.852, 
x Q » -0.0784, = 0.1165 and for F = -0.4 with the distur- 

bance application x(j) * x Q , x(j+l) =* x^ where j =» 2nIFHa, 
n = 5, m = 0,1,2, ..... 7. respectively. Here again the 
. transients between two steady state solutions are observed 
when the disturbance is applied at various die arete time 
instants. 

For V = 5,6, the results obtained are verified with 
exact simulation and the time plots showing the jump behaviour 
are not given here. 



mm 


(a) Jump between solutions x (40)=~0*0784, x (41 )= 0-1 1 65 


(b) Jump between solutions x(41)=-0*0784,x(42)=0 / l165 


W 


(c) Jump between solutions x (42 )=- 0-0784, x (43) = 0-11 65 


warn 


(cl) Jump between solutions x(4 
FIG-5-8 SOLUTION TO-EQN 15-ii 


x (44)= O H 65 


b*-0'95 , F* 0 






x ( K) 0-0 


10 


80 


- 1-0 


(e) Jump between solutions 



x(44M>0784,x (451=0-1165 
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She region ABCS! in Big* 5*3 on the otherband satis- 
fying 0-2 and 0-3 alone, sustains nonlinear oscillations. 

Ibis oan he verified by picking coefficient values £a,bj 
Inside this region and simulating the system given in (5*1) 
for initial values given in (5.6) and (5.7). This response 
is shown in Fig. 5*9* She nonlinear sustained oscillation is 
observed without any jump for all values of k with small ini- 
tial transient. Fig. 5.10 provides the time response of 
non! inear sustained oscillation for parameter values selected 
from the region where such oscillations are possible for V*4. 

As mentioned earlier the xange of Input amplitude F for 
which there exists a jump can be obtained as follows. For • 
specified filter coefficients [a»b] and for input frequency 
a = w/g, the value of F satisfying all the three Inequalities 

0- 1, 0-2 and 0-3 gives the required range. For example for 

1- ?, a - 1.95, b - -0.995, the F satisfying the 0-1, 0-2 
and 0-g conditions in sqn. (5.18) cure 

F » -0,95251 
F > -1.27105 
F < -0.76214 

respectively. 

Then the range of F satisfying all the above three 
inequalities is given by. 


—0.76214 ■ ^ F < —0.95251* 


(5.42) 
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Selecting F = -0.8 (which is within the range given in 
eqn, (5.42)) the time response is given in Fig. 5*11. Fig, 
5.11a shows the linear response with x(o) = 0.426 and 
x(l) = 0.840, Fig. 5.11h gives the nonlinear oscillation for 
x(o) = x(l) = 1,0 and Fig. 5.11c shows the jump between the 
solutions , 

Then for N =.4, a =1.95, b = -0.99, F satisfying the 
conditions given in eqn. (5.22) are 

F > - 0,5429 

F > 0.5462 

F < 0.4222 

from which the range of F satisfying all the above inequa- 
lities is given by 

-0.5429 < F < -0.4222 

The plot shown in Fig. 5.12 gives the time response of the 
system for the coefficient values given above and with 
F = -0.446. Here again the jump up behaviour between two 
steady state solutions is evident. 

Similarly for H = 5, a = 1.9,b = -0.95, the range of 
F is given by 

- 0.3116 < F < - 0.2746 

and for N =6, a =1.9, b = -0.95, the range of F is given by 

- 0.2032 < F < - 0.1889 
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Figs. 5.13 and 5*14 provide the time plot indicating the 
jump for F = -0.274 and F = -0*188 for N = 5 and N = 6 
respectively. 

5.5.2 Jump down phenomena; 

In the previous section a jump up phenomenon between 
two steady state solutions is observed in digital filters 
under forced situation. Jump down phenomenon is also 
possible under the conditions discussed below. 

In this oase the solution is initially obtained for 
x ( o ) = x Q and x(l) = x 1 which is periodic with unity mag- 
nitude (nonlinear oscillatory solution). When the condi- 
tions given in eqns. (5.6) and (5.7) are applied such that 
x(j) = x(o) and x(j+l) = x(l), where j is the discrete 
time instant given by j = 2nIT for n = 1,2,3, a jump 

down phenomena is observed without transient as shown in 
Fig. 5.15a for N = 3. On the otherhand when the conditions 
x(j) = x(o) and x(j+l) = x(l) are applied at time instants 
other than the above, the output sequence shoots up to 
nonlinear mode after a few cycles of transient behaviour. 

This is indicated in Figs. 5.15b to 5.15e for IT = 3 with 
disturbance application at various time instants as indicated 
therein. Similar behaviour is observed for other frequences 
namely N = 4,5 and 6. 
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It is of interest to refer to a similar although, 
less complete investigation due to Jsjristianss on [ 95.] concerning 
the nonlinear oscillations m digital filters. As one readily 
sees, the proposed technique provides information regarding the 
nonlinear oscillations only for positive values of the parameter 
’a’ . For negative values of ’a* no jump phenomenon of the 
hind described is observed. This is, because there exists 
no F to satisfy the general necessary conditions C-l to C-3 
for specified N. In [ 95 ] it has been demonstrated the 
existence of jumps for small positive values of F (for M = 3), 
where as in the proposed analysis the admissible value of F 
turns out to be negative for N = 3 as seen m Fig. 5.5. This 
is because of a difference in the system model m [ 95 ] 
introducing an additional phase shift in input forcing f un ction. 

In the next section the occurrence of super/sub harmonic 
oscillations m a forced second order digital filter* is dis- 
cussed, 

5.6 Super/ sub harmonic oscillations : 

In nonlinear systems, the steady state response to a 
harmonic signal may have a frequency which is multiple or 
submultiple of input frequency and they are m general known 
as super and sub harmonic components. This kind of phenomenon 
in a nonlinear discrete time systems has been studied in 
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detail m chapter three for polynomial type of nonlinearity. 

In continuous tine systems the super/sub harmonic oscillations 
have been discussed extensively in many text books and research 
papers. Por the present study the available methods can not 
be extended for studying the super/ sub harmonic oscillations 
m nonlinear digital filters due to discontinuous nature of 
the nonlinearity. Attempts, however unsuccessful, were made 
to investigate the super/sub harmonic oscillations m non- 
linear digital filter with different known techniques, namely 
the harmonic balance method and two variable expansion 
procedure* Then the possibility of occurrence of super/ sub 
harmonic oscillations are deducted based on exhaustive simu- 
lation of given system. Two kinds of sub harmonic oscilla- 
tions are observed and the first one is the type, which has 
been reported by Claasen and Eristiansson [96], the time 
response of such oscillations shows that the sub harmonic 
components die out for large tine, whereas the second type 
of sub harmonic oscillation is a steady state oscillation. 

Pigs. 5.16 and 5.17 give these oscillations. Por these 
investigations the system description is same as given in 
[96]. 

The existence of nonlinear behaviour such as limit 
cycles, jump phenomenon and super/sub harmonic oscillations 
are generally not desirable in designing digital filters 
[115]. It is important therefore to find methods that avoid 
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as much, as possible, the occurrence of nonlinear phenomena. 
Glaasen et al [96] and Willson [97] have proposed error 
feedback circuits to improve the behaviour of digital filters 
minimising the nonlinear effects. It has been shown by them 
that the proposed feedback circuit do not completely eliminate 
the occurrence of overflow phenomena, but with optimal choice 
of feedback coefficients the nonlinear effects are reduced to 
a minimum. 

5.7 Conclusion : 

A second order digital filter with overflow saturation 
nonlinearity is analysed under forced situation. Some 
conditions on the filter coefficients [a,b] are obtained for 
given input frequency interns of the input amplitude for 
which the filter exhibits nonlinear overflow oscillations, 
such as nonlinear sustained oscillations and jump phenomena. 
a simple graphical construction of the derived conditions 
gives the region m the a-b plane (for specified input ampli- 
tude) in which either sustained oscillation or jump behaviour 
is possible. The range for the input amplitude is obtained 
for specified input frequency and parameter [a,b] values 
for which the nonlinear behaviour is observed. The occurrence 
of subharmonic oscillations are studied through exhaustive 
search by simulating the system equation. The theoretical 
results obtained are verified by computer simulation. 



CHAPTER 6 


CONCLUS ION 

In this thesis a detailed analysis of a class of 
nonlinear difference equations arising generally in dis- 
crete time systems has been carried out. The nonlinear 
difference equations resulting m a discrete time system 
are generally of two kinds, 

i) Those described by polynomial type of nonlinearity and 
11 ) Those described by discontinuous type of nonlinearity . 

The second type of equations generally exist in the 
mathematical description of digital filters under nonlinear 
operation due to wordlength reduction that is necessary on 
account of finite wordlength registers. In this thesis these 
two kinds of nonlinear equations have been analysed without 
and with external excitation. Nonlinear phenomena like limit 
cycle oscillations under force free condition, jump phenomena 
and response characteristics with weak forcing and sub/ 
superharmonic solutions under strong input situation are the 
mam investigations. 

The difference equations described by polynomial non- 
lmeanties have been analysed adopting Discrete Multiple 
Scale Perturbational Technique. This technique has already 
been discussed m depth by Cole [16], Nayfeh [18] and others 
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for continuous time systems. This technique has been so far 
applied to study a class of ordinary and partial differential 
equations as well as a class of third order differential 
equations [21], Yen and Krouauer [20] have applied the two 
time variable expansion procedure to analyse three oscillators 
coupted by a weak nonlinearity, hronauer and Musa [25] have 
studied the possible sub/superharmonic oscillations in a class 
of second order differential equations using this technique 
under strong single input situation. Whereas Tiwari and 
Subramanian [22] have considered the two time scaling technique 
to investigate super/sub/ultrasubharmonic solutions m a class 
of nonlinear second order differential equations with strong 
multiple inputs. The method adapted in this thesis, namely 
the Discrete Multiple Time Perturbational Technique has been 
derived based on the properties of finite difference operators. 

The proposed method has been tested initially for linear 
difference equations with under damped, critically damped and 
unbounded situations and the analytical results are compared 
with simulation results. Dor free and weakly forced nonlinear 
systems, a discrete version of a Duffing oscillator has been 
considered and the well known response characteristics, loci 
of vertical and horizontal tangencies have been studied both 
with and without damping. The existence of limit cycle 
oscillation has been shown by considering a Van der Pol type 
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of second order difference equation and the stability analysis 
of such limit cycles has been carried out through a variational 
technique. In order to reduce the mathematical calculations 
involved m the proposed method, a modified definition for the 
two time discrete model has been provided and the results 
obtained by this modified technique have been compared with the 
other approaches. 

Strongly driven second order (weakly nonlinear) difference 
equations have also been considered and generation of sub/ 
superharmonics has been investigated. Both the discrete 
multiple time perturbational technique and the harmonic 
balance method have been adapted in the above investigation 
and the two results compared. 

Free second order digital filters with saturation over- 
flow nonlinearity have been considered and the existence of 
limit cycle oscillations as well as a monotonically decaying 
type of response for parameter values outside the stability 
triangle in the parameter plane have been demonstrated. 

Regions in the parameter space xtfhere the limit cycles of 
different periods are possible and locations both in the 
parameter plane and in the initial condition plane where the 
output sequence decreases monotonically have been mapped. In 
the analysis of limit cycle oscillations three methods have been 
proposed for locating the various pertinent regions m the 
parameter plane. The first of these is a variation of one 
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proposed by Ebert et al [76] whereas the other two arc based 
on a polynomial approximation to the saturation nonlinearity. 
The third method exhibits the use of discrete multiple time 
perturba „ional approach and is believed to be quite novel. 

Harmonically forced second order digital filters have 
also been considered and the existence of jump resonance, 
super/ subharmonic solutions have been studied with saturation 
overflow nonlinearity. Knowing the input frequency and the 
filter coefficients, conditions have been derived for the 
nonlinear oscillations to exist. A range of values for the 
input amplitude has been obtained for the existence of non- 
linear oscillations, such as jump phenomena and nonlinear 
sustained oscillations. The minimum perturbation or distur- 
bance (from a linear oscillatory solution) necessary to cause 

the above said nonlinear oscillations have also been obtained. 
The results obtained for jump up and jump down conditions have 
been verified by computer simulation. 

Scope fo* Further Work ; 

The problems listed below need further investigations 
m the study of free and forced oscillations m a class of 
second order nonlinear difference equations with polynomial 
as well as discontinuous type of nonlmearities . 
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1. The Discrete Multiple Time Perturbational Technique 
proposed m this thesis appears to he nonunique m the sense 
that the derivation is based on the properties of finite 

diff erer oe operators. It would be worthwhile investigating 
whether other techniques based on theory of finite difference 
operators are available for analysis of discrete time systems. 

An appropriate comparison can be made. 

2. It would be interesting to work out some general exi- 
stence conditions for presence of nonlinear phenomena like 
limit cycle oscillations, jump phenomena etc. 

3. The proposed multiple scale perturbational technique 

can be extended to study certain higher order as well as coupled 
nonlinear discrete time systems. 

4. The scheme 2, (introduced in section 2.8 of this thesis) 
namely the modification over the proposed method has been used 
to obtain more accurate solution for equations whose response 
is known to be bounded. It would be interesting to investigate 
m detail as to tho reasons for which the modified definition 
m eqn. (2.90) gives more accurate qualitative results. 

5. For force free second order digital filter a region m 
which a monotonic type of decaying response possible has been 
investigated. It would be interesting to investigate existence 
of other kinds of decaying responses for instance, oscillatory 
decay. 
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6. In forced second order digital filter with saturation 
nonlinearity the existence of subharmonic oscillations have 
been investigated by exhaustive simulation of the system 
equation. It would be useful to propose an analytical 
technique for investigating the above phenomena m a discrete 
system described by discontinuous type of nonlmeanties, 
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APPENDIX A 

PROPERTIES OP THE STABILITY TRIANGLE 

It is well knowxi that fur a linear second, order difference 
equation, 

x(k+2) + ax(k+l) + bx(k) = 0 (A.l) 

the stability of the response x(k) can be studied knowing the 
location of the system parameters 111 fke a ~k parameter 

plane (shown in Eig. A.l). When [a,b] values are inside the 
triangle AOG-, the system response is asymptotically stable 
(decaying to zero), when they are on the sides of the 
triangle ACG-, the solution x(k) is generally bounded (periodic) 
and for the parameter location outside the triangle aOG the 
response is always unstable (unbounded response). The tri- 
angle aCG is known as 'the stability triangle*. 

It is to be noted that the nature of the response of 
system m eqn. (A.l) depends on the location of [a,b] values 
in the parameter plane. The response x(k) for different 
location of system parameters are given below. 

(1) at A, 

x(k) = a + Bk, unbounded response without oscillation. 
x(k) = (A+Bk) (-1)^, unbounded oscillatory response. 


(2) at C, 
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(13) along the curve CE, 

cL 

x(k) = (a+BJc) (- ^-) , oscillatory decay. 

(li| ) region I, 

x(k) = (Vb)^ (a cos G^k +• B sm Q^k), oscillatory 
decay. 

where 9^ = cos~' I "(a/2/b) . 

(15) region II, 

x(k) = (Vb)^ (a cos Q^k + B sm G^_k), oscillatory 
decay. 

where 9^ = rc - cos - " 1 ' (a/2ft>). 

(16) region III, 

x(k) = a(a)^ + B(p)\ decays to zero solution 
a and j3 are positive, less than unity. 

(17) region 17, 

x(k) = + (-p)^, oscillatory decay. 

a and p are positive, less than unity. 

(18) region 7, 

x(k) = A(a)'^ + 3(-p)^, decays to zero solution 
a and p are positive, less than unity, a > p. 

(19) region 71, 

x(k) = a (- a)^ + B(p) ic , oscillatory decay, 
a and p are positive, less than unity, a > p. 

(20) region VII 

x(k) = A (a)' 11 + B(-p)^', unbounded response. 

a > 1.0, 


0 > p > 1.0. 
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(21) region VIII, 

1^» lr- 

x(k) = ii(— a) + B ( p ) , unbounded oscillatory, 

a > 1.0, 0 > p > 1.0. 

(22) region IX, 

x(k) = A ( a ) ^ + B((3) , unbounded response, 

a > 1.0, 0 > (3 > 1.0. 

(23) region X, 

x(k) = a (- a)^ + B(-p)^, unbounded oscillatory, 
a > 1.0, 0 > p > 1.0 

(24) region XI, 

x(k) = (fb)^ (a cos 9^k + B sin Q^k), unbounded 

response, b > 1.0 (the region XI is inside the 

2 

parabola a = 4b). 

where A and B are the constants to be determined from 
the given initial data. 
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APPENDIX B 

DISCRETE MODEL FOR A DIFFERENTIAL EQUATION 

The discrete model of a given differential equation is 
obtained by application of the central difference operator 6. 
The central difference operation is defined as follows. 

Let x(k) be a real function of integer variable k, the 

2 

first and second central difference operators 6 and 5 are 
defined as 

6x(k) = x(k+-i-) - x(k— i-) 

6 2 x(k)= x(k+l) - 2x(k) + x(k-l) (B.l) 

The two differential equations for which the discrete 
models are obtained are the Duffing' s equation with and without 
damping. These are indicated below as equations (B.2) and 
(B.3) respectively. Note that since there may exist a phase 
difference between x(t) and the impressed force (due to the 
presence of damping term)* for the Duffing equation with 
damping [5], which is other than 0 and % f the model is assumed 
to be of the particular form shown in eqn. (B.2). For situa- 
tion where there is no damping, model m equation (B.3) is 
adequate. 

x + pcx + ax + pyx"^ = p(F^ cosw t - sm w t) (B.2) 


• * Z 

x + ax + pyx = pF cos w t 


(B.3) 



285 


where ■ refers to a, c, y are constants and (i is small 

parameter. 

The discrete models for equations (B.2), (B.3) are now 
obtained by application of < equation (B.l) rs 

x(k+l) + X x(k) + x(k-l) + [ic [x(k+f) - x(k~ 0] + [iyx^(k) 

= [i[P^ cosw k - F 2 smw k] (B.4) 

and 

x(k+l) + X x(k) + x(k-l) + }iyx^(k) = [iF cos w k (B.5) 
where 

X = a - 2. 
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APPENDIX 0 

POLYNOMIAL APPROXIMATION TO A SATURATION NONLINEARITY 

a polynomial approximation for a given saturation 
nonlinearity is possible for a given restriction on the range 
of the dependent variable. The saturation nonlinearity and its 
mathematical representation are given below m Pig. 0-1 and 
m eqn. (C-l) respectively. 

f(cr)t 


-l.( 

/ 

Pig. 0-1 Saturation Nonlinearity 

J > 1.0 

| orf < 1.0 

cr < -1.0 (C.l) 

Since the nonlinear characteristic curve shown in Pig. 
0-1 is confined within first and third quadrants, the poly- 
nomial approximation takes the following general form. 
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f(cr) = a-j^a + a^cr 5 + a^cr 5 + .... (C-2) 

where are the coefficients to he determined so that the 
approximation given m eqn. (C-2) is as close as possible 
to the given saturation nonlinearity . 

The coefficients a are determined, by assuming a range 
D for a, the degree of the polynomial, considering as many 
points as possible m the range 0 to D and computing using 
the pseudo inverse approach, as follows. 

Assume the degree of the polynomial to be used as m 
(m odd) and the range of a be D, that is 

I a I < D 

Then consider n equal points in the known range 0 to 

D, namely a- L , a 2 , cr^ c n . With this construction, the 

following set of equations are obtained from the eqn. (C-2). 

f (cr^) = a-^ + a^o^ + + + a m <3^ 

f(ff 2 ) = a x 0 2 + a 3°2 + a 5 0 2 + ***•• + 

f(oj) = a 1 o 3 + + a^ + + a m o^ 
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Eqn. (C-5) can be rewritten as follows in matrix form : 

a x = (C-4) 

where, 

n i n 

A is a known rectangular matrix of order nx(=-j=-) and 
is given by 



and y is a known n vector 
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f < a l> 

f (<J 2 ) 
f(cr 3 ) 


f <V 


m 

Then multiplying (0-4) hy A (transpose of A), we have 


T T 

a A x = A £ 

x = (A T A) -1 A T y (0-5) 

eqn. (C-5) gives the values of the coefficients of the 
polynomial in eqn. (Cg). 

For example, with m = 7, D = 3.0, n = 12, the following 
values for a^ are obtained ; 

= 1.16253 
a 3 = -0.2850 
a^ = 0.0364-3 
a rj = -0.00167. 

Thus the polynomial approximation is 

f ( cr) = 1.16253 o - 0.285 + 0.03643 o 5 + 0.00167 cr 7 


with 


| cr f < 5.0 • 
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